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TECHNICAL PAPER 


REVIEW OF THE PROBABILISTIC FAILURE ANALYSIS METHODOLOGY AND 
OTHER PROBABILISTIC APPROACHES FOR APPLICATION IN 
AEROSPACE STRUCTURAL DESIGN 


I. INTRODUCTION 


The purpose of this report is to document the work the Marshall Space Flight Center 
(MSFC) probabilistic analysis team has been involved in during 1991-1992. The report introduction 
is divided into three sections: (A) Background Information; (B) MSFC Task, Objectives, and 
Operating Plan; and (C) Project Scope. 


A. Background Information 

For a number of years, NASA has been designing spacecraft and structures using the safety 
factor/deterministic approach to design. Although the method is well established and generally 
accepted by the aerospace community, it does not provide a quantitative means for identifying the 
failure risk associated with a given design or flight. In 1985, the Jet Propulsion Laboratory (JPL) 
Engine Certification Project was initiated to develop an improved methodology for quantitatively 
evaluating and establishing flight readiness (certification). The program was titled Probabilistic 
Failure Analysis (PFA), and it is funded through codes Q and M. The Challenger accident (1986) 
impacted the development of this technology in several ways. First, Professor Richard Feynman, 
during the failure investigation, determined that the design engineers believed the risk of engine 
failure was about 1 in 200; higher management understood the risk to be 1 in 1 00,000. 1 Later analy- 
ses supported the engineers estimate. Somewhere between design and certification for design, the 
real information was being lost. NASA began to examine different approaches for identifying risk, 
perhaps at the design level. Risk could then be elevated to certification in a quantitative procedure, 
rather than qualitatively implied from safety factor. 

The second influence of the Challenger accident was that a significant amount of money was 
appropriated for the advancement of safety/risk technology. In addition to the JPL Engine 
Certification Project, the Probabilistic Structural Analysis Methodology (PSAM) Project, directed by 
NASA Lewis Research Center (LeRC), received large sums of funding. The primary analysis tool 
that resulted from this project was NESSUS, a probabilistic finite element software code developed 
by Southwest Research Institute. The PSAM program was structured as follows: (1) to understand 
how probabilistic analysis fit into the engineering process at the component level; (2) to advance 
probabilistic analysis into the system level; and (3) to develop the understanding of how it fits into 
the certification process. In other words, the PSAM Project approached the safety/risk technology 
problem from the bottom up, while the JPL Engine Certification Project approach was from the top 
down. Although the philosophies vary, both methodologies handle problem uncertainty solely from a 
probabilistic format, not a safety factor approach, or any combination with safety factor. 
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In 1988, an effort was initiated with MSFC to begin “technology transfer” of the JPL prob- 
abilistic methods. The team was led by Henry Stinson of the Propulsion Lab. From 1988 through 
1990 little progress was made. The PFA technology, as defined, documented, and presented, was 
extremely difficult to penetrate and understand. In January of 1991, 6 years after the JPL Engine 
Project was initiated, no one (except the JPL team) was applying the PFA technology to real prob- 
lems. As a result, the PFA program was canceled by the Shuttle Program Office (Bob Crippen). 
Later, the PFA program was reinstated with limited funding through FY92, provided MSFC would 
take on an active roll. In March of 1991, a new Marshall team was established. The results and 
findings of this MSFC team’s efforts (May of 1991 to December 1992) are documented herein. 


B. MSFC Task, Objectives, and Operating Plan 

The MSFC probabilistic analysis team consists of the following members: Dr. John 
Townsend, Mario Reinfurth, Rene Ortega, Charlie Meyers, and Jeff Peck, Structures and Dynamics, 
and Bob Weinstock, Vitro. A number of others participated during the investigation. Namely, Greg 
Swanson from Stress, Brian Goode from Thermal Division, and Steve Gentze from Materials. The 
JPL team participants are Dr. Nick Moore, Dr. Don Ebbeler, Ms. Laura Newlin, and Dr. Dharshan 
Sutharshana. 

A program operating plan was jointly written by the MSFC and JPL teams and approved by 
Headquarters codes Q and M. The plan is documented in reference 2. The overall program objective 
was to formally assess the merits and application of the JPL team’s PFA methodology in the design 
process. No formal assessment of the PFA methods for use in certification was made. The program 
goals were: (1) to develop an indepth understanding of the PFA methodology, (2) to assess the 
utility of the methodology and supporting software tools, and (3) to develop a plan for the technology 
transfer of the methods from JPL to MSFC. From the outset, the MSFC team decided not to use the 
PFA methods as a “black box” tool into which information is input and out pops an answer. Also, 
the team chose problem application as the basic approach for achieving these goals. To climb the 
learning curve as quickly as possible, a simple cantilever beam example problem was worked using 
the PFA software and approach. This simple beam example provided an excellent tool for communi- 
cation between JPL and MSFC, as well as to others interested in understanding the basic PFA 
methods. The other problem worked of interest to the Marshall team was the high pressure oxygen 
turbopump (HPOTP) first-stage turbine disk seal rib cracking problem. A third problem identified as 
the carrier vehicle support structure — part of the Aeroassist Flight Experiment (AFE) — was 
removed from the task, since AFE was canceled in early 1992. Documentation review was also a 
major part of the MSFC task team approach for penetrating the PFA methodology. 

C. Project Scope 

The scope of this report is divided into several sections. Section II documents the details of 
the cantilever beam example. In particular, this section identifies: (1) the steps involved in working 
this type of problem using a probabilistic format and (2) findings in terms of the PFA approach com- 
pared with a more classical probabilistic analysis. Section III presents specifics of the PFA methods 
and, in particular, how these methods fit into a probabilistic framework. Also documented in section 
III is an issues and concerns review and a class-of-problems list for which the PFA approach has 
application. Section IV discusses the methods and software developed by the JPL team to model the 
material uncertainty characteristic of S-N fatigue data. The classical approach for handling fatigue 
data uncertainty is also examined. Section V provides the reader with insight into the application of 
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response surface methods; detailed examples are presented. Also, the role of response surface 
methods in the PFA process is reviewed. A detailed discussion of the HPOTP problem example is 
given in section VI. Other methods, and perhaps better approaches at using probabilistic methods in 
engineering, are presented in section VII. Of particular interest is a review of the finite element 
probabilistic engineering analysis tool known as NESSUS. Section VIII gives a short discussion 
with recommendations on the development of a NASA probabilistic-based design code. The 
document is concluded in section IX with a report summary and recommendations for application of 
probabilistic methods within NASA. 


II. CANTILEVER BEAM EXAMPLE 


A simple cantilever beam example problem was chosen as the educational and communica- 
tion tool of the MSFC team. The beam problem helped the team to rapidly climb the learning curve of 
probabilistic analysis concepts and, in particular, the PFA methods. Also, it provided the MSFC 
team with a communication tool for demonstrating these concepts and findings to others. 


A. The Steps 

Figure 1 gives a sketch of the cantilever beam problem along with the steps used to work the 
problem using a probabilistic format. The first step is to develop a “closed-form” equation that 
defines the failure parameter as a function of the model parameters. In this case, beam tip displace- 
ment is defined as the failure parameter with input parameters, termed drivers, of load (P), beam 
length (L), material modulus of elasticity (E), and cross-sectional moment of inertia ( I z ). To develop 
a closed-form equation of this nature is not trivial, it requires engineering knowledge of problem 
physics, and in many cases simplifying assumptions. For example, to define the equation shown in 
figure 1 requires that one assume linear and homogeneous material behavior, small displacements, 
and no shear deformation. 

The second step in the process is to establish the model parameters or input probability dis- 
tributions. This part of the procedure is probably the most difficult for engineers to complete. Three 
methods for obtaining the distributions are: (1) test data, (2) design requirements, and (3) engineer- 
ing judgment. The preferred model is the one based on test data. However, in many cases, test data 
may not be available or economically obtainable. If the engineer is going to use a distribution based 
on his judgment, then it is a good idea to have a design tool (software) available for establishing 
parameter distribution sensitivities. 

Once the math model is determined in closed-form, and the driver distributions are chosen, 
the third step in the process can be completed. Namely, the engineer can use Monte Carlo simulation 
to determine the probability distribution of the failure parameter/output. Several different computa- 
tional methods are available for estimating the reliability of a design, and simulation is perhaps the 
easiest understood. Also, simulation is the basis for all of the PFA methods. The procedure is as 
follows. A random number generator is used to pick a single value from each of the input parameter 
distributions. These numbers are then used in the closed-form equation to give a single deterministic 
value for the output parameter. This same loop is repeated thousands of times until enough informa- 
tion is available to either determine the probability of nonexceedence (estimated from counting) or to 
fit a distribution to the output parameter data. 
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Figure 1. Cantilever beam example. 
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The final step addresses the question of how to use the data. Two possible avenues are 
available. The distribution data can be used (1) to estimate structural reliability, or (2) to understand 
output/input parameter sensitivities. First of all, there are two ways to estimate reliability from data. 
The simplest and most straightforward method is to count the number of times the output parameter 
exceeds the failure criteria. For example, in this cantilever beam problem, the failure criteria was a 
tip displacement of 1 inch. The procedure is to run the simulation X number of times (say, 5,000) and 
count the number of times the answer exceeded 1 inch (say, 5). This direct method gives an estimate 
of reliability of l-(5/5,000) or 0.999. The counting method does not require curve fitting or 
extrapolation; but, typically, it does require a large number of computer runs. 

The second approach for determining the reliability estimate is through curve fitting to either 
the entire output distribution, or to the tail region of the distribution. Of course, this procedure does 
not require as many computer runs, but it does use the curve fit for extrapolating data. Many 
engineers consider data extrapolation to the reliability estimate an absolute wrong. The problem is 
complicated by the fact that math models and parameter distributions used to perform the analytical 
estimate of reliability are in error. Furthermore, the engineer must be prepared to make confidence 
statements with respect to the reliability estimate. 

By definition, true reliability is demonstrated, not simply estimated from an engineering 
analysis. Until failure rate data bases are available from experience, probabilistic methods can best 
be utilized as a design tool to help identify the sensitivities of problem parameters. 


B. The Findings 

The input drivers in the problem are not shown graphically, but the distributions used, their 
mean values, and the corresponding coefficients of variations (COV’s) are referenced below: 

INPUT DATA 


DRIVERS 

BEAM LENGTH (L), inches 

MOM OF INERTIA (IZ), inches 

APPLIED LOAD (P), kips 
APPLIED LOAD (P), kips 

YOUNG’S MODULUS (E), lb/in 2 


PARAMETERS 

UNIFORM (34.200, 37.800) 

UNIFORM (3.665, 5.470) 

NORMAL: (1,000.0, 25.0)-Figures 2 and 3 
LOGNORMAL: (1,000.0, 25.0)-Figures 4 and 5 

NORMAL: (0.30E+08, 0.60E+06) 


Results of the cantilever beam example problem are shown graphically in figures 2 through 5. 
Notice in figure 2, the probability density function (PDF) for the output parameter, beam tip 
displacement, is plotted by the top graph; and the cumulative distribution function (CDF) which is 
the area under the PDF curve (i.e., the integral of the PDF function) is given by the bottom graph. 
Three standard probability curves are used to fit the simulation data (5,000 values); the normal, the 
Weibull, and the lognormal. The sample mean of tip displacement is 0.116 inch, with a standard 
deviation of 0.017 inch and a COV of 14.8 percent. Beam tip displacement is plotted on the graph 
abscissa, and the PDF value is the ordinate. The data points define a histogram grouping of 5,000 
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Figure 3. Cumulative distribution function curve fits of cantilever beam displacement data. 
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BEAM TIP DISPLACEMENT, 1/X , INCHES 


Figure 5. Left-hand tail of Burr distribution versus simulation data — cantilever beam problem. 

simulation values. Twelve bins are used with a bin width of 0.01 in. Also note, the histogram data 
have been normalized so that the area under the curve is one. Keep in mind that a different choice of 
bin size and number will result in a slightly different representation of the data values. The PDF 
curve fits to the simulation values are also approximations. Of course, the truest representation of 
the simulation data is given by the data values (i.e., no curve fit). 

Based on visual comparisons, the lognormal curve fit provides the best analytical data 
approximation. This finding is easily verified using statistical curve fitting routines and “goodness- 
of-fit” tests. The MSFC team used a program developed by Professor Paul Wirsching to determine 
which curve provides the best data fit. 3 The program uses a W-statistic based on a form similar to 
the Cramer-Von Mises statistic. If engineers are going to use probabilistic methods to analyze 
structures, then it is imperative that the probabilistic design tools include methods for determining 
which curve or curves best represent the input drivers and output data. The current PFA software 
package does not included “goodness-of-fit” tools. 

The cumulative distribution curve, known as the CDF or Ogive, is the characteristic curve of 
probabilistic analyses. Figure 3 shows an expanded view of the right-hand tail of the CDF curve 
given in figure 2. Probability predictions based on normal, Weibull, and lognormal curve fits are 
shown, along with a histogram reduction of the data (identified by the CDF data points). Also shown 
in the graph is a probability curve based on counting. For example, each point on the curve defines 
the quantity 1 minus the ratio of the number of exceedences of that value to the total number of 
simulations. The counting procedure, of course, is the preferred method of estimating reliability; but 
typically, it requires a large number of computer runs. 
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Since in a well-designed system the probability of failure is low, the engineer is concerned 
usually with the extreme values of the input drivers (design variables) and tail ends of the output 
parameter distributions. For the cantilever beam example given, the corresponding probability 
number to a tip displacement of 0.16 inch is highly dependent on which type distribution is used to 
represent the simulation data. A Weibull curve fit over estimates the reliability at a value of 
0.9999999, and the lognormal curve fit under estimates the number at a value of 0.985. The counting 
method gives a value of 0.9986. The engineer must be concerned with not only developing a good 
model, but in understanding how the design variables influence the reliability estimate. Also, if the 
engineer attempts to curve fit the output distribution, then he must determine if the final reliability 
estimate is conservative or unconservative, and to what degree. Until reliability engineering 
experience bases are developed for a particular material and design, it is highly recommended that 
the counting methods always be used to estimate reliability. Data extrapolation is not recommended, 
and interpolation methods must be used within their range of accuracy. In any curve fitting procedure, 
convergence limits must be defined. The bottom graph in figure 3 shows the effects on probability 
estimates of using a Weibull curve fit of 5 simulations versus 5,000 simulations. 

The results of a second cantilever beam example problem are graphed in figures 4 and 5. This 
example shows some of the unique features of the PFA methods. The PFA approach uses the same 
basic steps in the design process identified earlier. The primary difference in the PFA approach is 
how the engineer is directed to handle the data. Figure 4 gives the PDF and CDF curves for the 
beam parameter of inverse tip displacement. To use the PFA method requires that the engineer 
transform the right-hand tail exceedence problem into a left-hand tail distribution. The transforma- 
tion process is necessary because the Weibull distribution curve, a left-hand distribution curve, 
forms the basis of the PFA methods. It does, however, complicate the engineering procedure and 
understanding of the problem. The engineer is concerned with tip displacement and its characteristic 
distribution, not with the inverse. Also, a distribution and its inverse and their means and standard 
deviations may not map the same. What is normal in one case can be lognormal for its inverse. 

The purpose of the transformation is to curve fit the simulation data to a Weibull-based dis- 
tribution known as the Burr distribution. Typically, a Weibull distribution is a two parameter distri- 
bution. The two parameters arc known as the characteristic life parameter, lambda, and the shape 
parameter, eta. The PFA method fits a two parameter beta distribution to the parameter lambda. 
This process maps the data into a three parameter distribution known as the Burr. Both the 
lognormal and Burr curve fit distributions of the PDF and CDF are plotted versus the normalized 
histogram data in figure 4. Notice, the Burr curve fit is only accurate in the extreme left-hand tail 
region of the curve. Figure 5 is an expanded view of figure 4. By tuning three parameters, the Bun- 
curve fit is matched closely to the simulation data to fit the left hand tail of the output distribution. 

The Burr curve fit is a very complex tool to use to model data. As it turns out, the purpose of 
the Bun curve fit is twofold: (1) it provides a good curve fit of the simulation data for extrapolation, 
and (2) it provides a math model that handles Bayesian updating in closed-form. The MSFC team 
believes both of these reasons have little application in practical design. First, reliability estimates 
based on extrapolation are not recommended. If curve fitting data and extrapolation are necessary, a 
simpler method is to use a three- or four-parameter least-squares curve fit to obtain the desired 
estimate. Second, Bayesian updating can only significantly change a reliability estimate (a prior 
estimate) if failure occurs. A “damaged” structure usually requires redesign. If a structure has failed 
and has been redesigned (typically, a costly process), then the logical and best approach to 
estimating risk of the new design is to develop new math models for use in updating the reliability 
analysis. Hence, Bayesian updating in this context has limited practical application. 
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in. THE PFA METHODS 


The JPL PFA methods and approaches to solving structural design problems are discussed in 
this section. Topics include (1) the probabilistic framework into which the PFA methods are 
structured, (2) a summary of the issues and concerns of the MSFC team with application of the PFA 
methods, and (3) the class of design problems the PFA methods are best suited for. 


A. Probabilistic Framework 

The block diagram shown in figure 6 identifies the characteristic steps of the PFA engineering 
process. The first step in working a probabilistic format, or for that matter any standard format 
problem, is to develop an engineering model of the structure or problem. The finite-element method 
is one of the primary analysis tools used by engineers for developing structural models. A typical 
finite-element model can have as many as 5,000 degrees of freedom (DOF). Some models are so 
detailed and complex, that 50,000 DOF is normal. Of course, this type of problem is more difficult to 
work than the closed-form problems such as the one defined by the cantilever beam in section n. 

There are two approaches for incorporating the finite-element model into a probabilistic 
framework: (1) one can simply use the entire model and perform Monte Carlo simulation, or (2) one 
can use the finite-element model to develop a response surface equation to use with the probabilistic 
analysis. The response surface approach is the only practical method for large models. The basic 
concept of response surface is to develop a single equation that can be used in place of the 
thousands of finite element equations. The engineer must make a lot of assumptions if he is to use 
this “optimizing” approach. The procedure is very similar to the design of experiments method for 
developing an empirical equation. First, the engineer must decide on the output parameter of interest 
and its input drivers. He must then write an equation in functional form; i.e., using a first-, second-, 
or third-order equation with coupling terms. The unknowns in the equation are the coefficients of the 
input driver variables. After the empirical form of the response surface equation is assumed, the 
finite-element problem is run a sufficient number of times to properly determine the coefficients. 
Further details of this method along with an example problem are given in section V. 

The response surface method provides a practical means for incorporating complex finite- 
element models into a probabilistic format or framework. The engineer must keep in mind, however, 
that a response surface is only an approximation to a finite-element model solution, which is also an 
approximation. It is important, therefore, to understand all the assumptions made in developing 
these two approximations, and the effects of these assumptions on the output parameter of interest. 
Once the response surface is defined (i.e., the closed-form equation is developed), the problem then 
reverts back to the same basic procedure used previously to work the cantilever beam example. 

A number of computational schemes are available for working probabilistic engineering type 
problems. The PFA methods developed by the JPL team are restricted to Monte Carlo simulation. 
Once the closed-form relationship is defined between the input and output variables, the procedure 
to use the Monte Carlo simulation is easily applied. Also, there are a number of software simulation 
tools on the market that can be purchased at a minor cost. Before simulation can be run, the input 
parameters, or life drivers, must be characterized in terms of a distribution. As mentioned in the 
beam problem, this step is probably one of the most difficult for the engineer to complete. In many 
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Figure 6. The JPL methods — probabilistic framework. 
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cases, no test data are available, and the engineer must make estimates on the type of distribution 
to utilize. It is important that the engineer understand the effects of his distribution assumptions in 
terms of sensitivity of input drivers on the output parameter. A good software tool is invaluable for 
these types of sensitivity studies. The PFA software is a FORTRAN coded package that is problem 
specific. The source code must be changed by the user to work a variety of problems. In this sense, 
the PFA software is not a good tool for the practicing engineer, but a research tool for use in 
development. 

The last major step of the engineering process is to generate simulated response data using 
Monte Carlo simulation of the response surface equation with the input life driver characterizations. 
The output parameter can be any parameter that correlates with life; i.e., stress or strain, or an event 
consequent parameter, such as burst pressure. The engineering process defined by these steps is 
not new technology. The key to making the probabilistic engineering process practical is being able 
to link it with established software tools such as finite elements. The response surface methods may 
provide the necessary link. The PFA methods make use of response surface technology, but the 
documentation is very limited and unclear. The assumptions that define a response surface must be 
well understood before application can be made. 

The main thrust of the PFA methodology is directed not to the engineering process, but 
toward the application of the output parameter data once they are generated. Three types of failures 
have been addressed; (1) event consequent failures (exceedence problems), (2) low- and high-cycle 
fatigue failures, and (3) fracture mechanics failures. Primarily, the MSFC team has reviewed the 
PFA low- and high-cycle fatigue studies. Fatigue data, in general, define a large scatter or 
uncertainty. Coefficients of variation (COV's) are typical at 50 percent. The method PFA uses to 
characterize fatigue data is addressed in the JPL report under the title of Materials Characteriza- 
tion. 4 It is a very difficult section to penetrate and understand. To help the MSFC team understand 
the process of characterizing the uncertainty in fatigue data, standard texts on the subject were 
referenced. 5 16 

Fatigue data are usually characterized by stress or strain life curves known as S-N curves. 
The standard practice is to place an uncertainty (distribution) either on (1) the life value ( N ) for a 
given stress or strain, (2) the stress or strain value (5) for a given life, or the slope of the S-N curve 
( m ). The distribution curve used to model the data should be defined by the data. In all cases, the 
PFA method uses a Weibull curve to represent the data. Once again, the reason for this fixed 
distribution choice is so that a closed-form equation can be developed for use with Bayesian 
updating. Also, the PFA material characterization method and software are structured with many 
features with limited utility for most problems. The standard probabilistic S-N concepts are not 
difficult to understand and apply. Section IV of this report describes the basic approach. 

The JPL team has done a considerable amount of work in developing the concept of 
bootstrapping for engineering application. The MSFC team has not reviewed their work in depth, but 
the bootstrap method may have application for characterizing fatigue data. The PFA fracture 
mechanics application is another area the MSFC team has not examined. A cursory comment is that 
the PFA fracture mechanics methods are structured to link with Bayesian statistical methods. Once 
again, the process is difficult to understand and apply. 

In conclusion, the MSFC team believes that the best way to introduce probabilistic methods 
into design engineering practice is through a well-established and defined process. The JPL team 
has failed to identify specific details of the engineering design process and how the PFA probabilistic 
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analysis methods apply. In particular, the technique by which an engineer integrates complex 
models, such as finite elements, into a probabilistic framework is very unclear. The PFA methods are 
structured to work primarily with the tail end of the engineering process; i.e., predicting a reliability 
number that can be later updated using Bayesian techniques. 


B. Issues and Concerns 

In reviewing the PFA methods, the MSFC team has a number of issues and concerns with 
these methods for their application to engineering design problems. Most of these concerns have 
been voiced earlier in the report. A summary list is provided below: 

1. Implementation: Finite Elements/Other Models . No specific guidelines on how to inte- 
grate complex failure models into a probabilistic framework have been documented, although 
response surface methods are recommended. If engineers are going to make use of reliability 
methods, then it is essential that the probabilistic analysis tools can be linked to existing analysis 
tools, such as finite elements. If the basic link is to be a response surface approach (design of 
experiments), then the engineer must understand what assumptions are fundamental to the problem. 

2. Life Driver Distribution Characterization . What distribution is the engineer to use for 
characterizing a specific parameter, such as the material modulus of elasticity? Curve-fitting routines 
and best-fit equations are not provided with the PFA methods. If test data are available, then the 
preferred method is to model the data with a best-fit distribution. However, in new designs, many of 
the data parameters are not test derived. Which distributions are recommended for these parameters 
and what are their advantages and limitations? Also, computer software is needed for checking out 
which distributions and parameters are sensitive to design. Basically, the PFA methods provide no 
specific guidelines or tools that help the engineer answer these questions. 

3. Probabilistic Design Computational Methods . The PFA methods are currently restricted 
to a Monte Carlo simulation computational scheme. The MSFC team agrees that simulation is the 
fundamental probabilistic analysis tool. However, it is also recognized that simulation is very time 
consuming and expensive for complex models. If PFA is to have any kind of application to practical 
engineering problems, then it must be structured with other computational methods that will speed 
up the analysis solution process. A number of other methods may have application. Also, “efficient” 
Monte Carlo simulation shows great promise. 

4. Computer Software . The PFA computer software is a FORTRAN coded “stand-alone” 
software package that is structured in modular fashion. The modular design gives it flexibility to 
work a variety of problems; but it is problem specific, and, therefore, requires reprogramming for each 
problem. This feature makes the PFA software package a good research tool, but a poor design tool. 
Essentially, the PFA programs should be restructured to make them user friendly. For example, the 
software should be structured to give the designer a menu of distribution choices for a parameter, 
and plots of these values. 

5. PFA Documentation . The reports that the JPL team have documented on the PFA 
methods are extremely difficult to penetrate and understand. The reports were written for those 
individuals that have a detailed understanding of probability theory and Bayesian statistics. They 
have not been written for use or application by the average engineer. This is a fundamental error of 
the JPL teams approach; i.e., not to adequately explain within their documents the details of how to 
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work problems using the PFA methods. A better approach would have been to show basic funda- 
mentals and to illustrate them with an example, so that an engineer not versed in these techniques 
could penetrate the theory and work a simple example to validate it. Documentation was structured 
for publication in technical journals, not for working engineers to understand. 

6. Statistical Methodology 

a. Limited to “Left-Hand” Tail Problems . The standard JPL methodology handles the 
case of simulating the left-hand tail of a Weibull distribution. However, many problems are not 
concerned with “early” failure, but with exceeding a fixed value. The focus is then on the right-hand 
rather than left-hand tail of the failure density. This is the case, for example, in the cantilever beam 
example where failure is defined as deflection exceeding Xo, where Xo is given. The current JPL 
approach is to “shift” the right-hand tail to a left-hand problem by considering the set {1/X}, where 
{X} is the set of simulated values. The apparent problem with this approach is that there seems to 
be no justification why the set { 1/X} follows a Weibull distribution. The mean and standard deviation 
values characterizing l/{ 1/X} and {X} do not agree; in fact the differences increase as the number of 
simulations increase. The PFA methods must be restructured to handle this set of applications (i.e., 
right-hand tail applications). 

b. Restricted to Burr Distributions . After obtaining a set of simulated failures, the 
statistical methods and FORTRAN programs used to determine reliability are restricted to a 
Weibull distribution where the parameter lambda ranges over a gamma distribution (known as the 
Burr distribution). The Burr distribution is just a method for fitting the very extreme left-hand tail 
region of the distribution using three parameters. Should not the choice of model distribution be 
governed by the simulated data? Of course, this approach would interfere with the “back-end” 
Bayesian statistics. However, the MSFC team believes the best approach to using probabilistic 
methods as a design tool is to not use a Bayesian logic. It makes more sense to model the data with 
the best distribution fit possible, especially for the purposes of data base development. 

c. Precision. Confidence, and Assurance Statements Confusing . A confidence level must 
be established for reliability estimations. After applying the JPL methodology to the cantilever beam 
problem and investigating a number of example problems, it appears that the JPL analysis tools are 
structured to handle confidence interval (or assurance as they call it) one specific way; i.e., to 
characterize confidence on the variable lambda as it ranges over a gamma distribution. A more 
general approach might be to put a confidence interval on the final simulated tail distribution using a 
nonparametric technique. 

d. Structured to Handle Bavesian Updating (Limited Application! . As the MSFC team 
understands Bayes updating laws, the probabilistic failure model results characterize the prior 
distribution, and the test data or flight data are used to update this distribution. In most cases, if a 
structure fails during qualification testing or flight, it is redesigned. A redesigned structure will be 
characterized by a new failure model, and the Bayesian statistics will not be applicable. For those 
cases where a failure does not occur, and Bayesian updating is performed, a less conservative 
reliability estimate will result — which may be applicable to life limit extension. From our viewpoint, 
Bayesian updating is not directly applicable to design problems, or problems that are structured to 
statistically characterize the likelihood of a specific failure mode. In those cases where we revisit the 
probabilistic failure model, the new information would be used to update or modify the driver 
distributions. The bottom line is that the MSFC team believes Bayesian updating is a very small and 
minor part of the whole engineering design process. 
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7. Materials Characterization (S-N Curves) . The methodology write-up for “materials 
characterization” needs to be clarified so that the underlying ideas are more easily understood. Also, 
over the past 50 years, numerous documents have been written on fatigue. And there are several 
methods for characterizing fatigue data uncertainty. Exactly, what is different about the PFA 
approach is difficult to understand. The bootstrapping methods are very interesting and should be 
better documented with examples. 

8. Model Accuracy Factor . The PFA methods use a modeling error factor that is uniformly 
distributed as an input or driver variable. What is the justification of assuming both conservative and 
nonconservative values; i.e., above and below 1? Would it not be more reasonable to use a constant 
value from the most conservative end of the distribution? In addition, it appears to us, that a better 
method for including modeling error and other inaccuracies is to incorporate them with the confidence 
regions of reliability at the end of the process. 


C. Design Application/Class of Problems 

The PFA methods are best structured to handle problems where the life variable can be 
expressed by a closed-form equation set. Also, if reliability predictions are required and the life 
variable defines an “exceedence type problem,” then the life variable must be transformed into a 
left-hand tail distribution. With respect to these two findings, a number of questions remain to be 
addressed; (1) What are the procedures to follow to determine a closed-form equation set from a 
complex finite element model (response surface approach is one alternative, are there others)? 
(2) What are the transformations that map into Weibull space for working exceedence type 
problems? (3) What are the underpinning assumptions in the PFA process? (4) What are the 
convergence criteria? 

Application of the PFA methods to sensitivity studies on driver variables offers the greatest 
advantage of the probabilistic approach to design. In high performance systems, such as rocket 
engines, it is extremely important that an engineer understand in detail the influence of driver 
uncertainties on the design. However, the PFA methods are not structured to automatically calculate 
the sensitivity values, a trait of a good design tool. The PFA process has been designed for 
application to specific fatigue and fracture type problems. Several specific examples have been 
completed by the JPL team. The engineer must keep in mind that any predicted reliability number is 
only an estimate, based on problem assumptions. Also, the accuracy of the reliability number is very 
difficult to determine. In no case should these analytical values be accepted, or used for 
certification. True reliability or risk must be based on data. 


IV. MATERIAL CHARACTERIZATION— JPL FATIGUE LIFE PROGRAM 
MATGRM/MATCHR, AN EVALUATION 


A. Stress/Life Statistical Analysis 

The purpose of this section is to describe the PFA fatigue model and how it can be duplicated 
using a much simpler method. The fundamental steps of least-squares curve fitting will be sum- 
marized, followed by an example. Additional background information is given in references 5 and 16. 
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The PFA statistical fatigue life model is basically a least-squares curve fit with a Weibull 
distribution on life rather than the commonly used lognormal distribution. 

The basic steps taken for a fatigue life model are as follows: 


(1) Collect raw data; stress, S and corresponding cycle life, N. Some data on Waspaloy, 
given below in table 1, will be used for this example. 


Table 1. S-N data for Waspaloy. 


Stress (lb/in 2 ) 

Cycles (jV) 

log ( S ) 

log (A0 

110,000 

39,600 

11.61 

10.59 

120,000 

11,349 

11.70 

9.34 

120,000 

3,749 

11.70 

8.23 

130,000 

4,163 

11.78 

8.33 

130,000 

3,824 

11.78 

8.25 

140,000 

4,743 

11.85 

8.46 

160,000 

1,019 

11.98 

6.93 

160,000 

677 

11.98 

6.52 

160,000 

636 

11.98 

6.46 


(2) Assume a relationship of the form. 


where. 


N = AS m , 


N = cycles to failure 
S = stress amplitude 
m = fatigue strength exponent 
A = fatigue strength coefficient 

which may be written in log form as 

ln(A0 = ln(/4)+mln(S) , 

which has the form of a straight line in log space, 

y = a+bx , 

where. 


y = log N 
x = log S 
b = m 
a = log A. 
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(3) a and b are estimated from the data using the least-squares method. This is a well- 
known approach described in most statistics texts. The following parameters are 
calculated from stress-cycle life data: 

S«=i(.x r x) 2 , 

1=1 

1=1 

n 

S v =Z(* r * ){y r y) , 

where. 


x, = ln(5,), log stress value 
x = ln(5), mean log stress 
= In (NJ, log cyclic life 


y = In(jV), mean log cyclic life. 

The following values are calculated from the data given in this example, 

5^ = 0.16035, 

5 W = 14.521 , 

Sxy = -1.419668. 

Least-squares curve fit constants are calculated directly as, 

5„ 

mean value slope, b = -s— = 1 12.7372 , 

^XX 

mean value intercept, a = y-bx = -8.85338 . 

(4) The least-squares fitted curve can now be compared with individual data points to 
determine whether or not the curve fit is adequate over the stress range of interest. A 
plot of the data with curve fit in log space is shown in figure 7. 
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S/N Data with Curve Fit 



6 7 8 9 10 11 


log(N) 

Figure 7. Waspaloy S-N data plotted with curve fit. 

(5) Next, determine estimates of the mean and variance of the cycle life for a given stress 
level. 


mean value cycle life, y = a+bx , 


2 Syy-iSfy) /£.„ 

variance, a A = — == 

n- 2 


The mean and variance are statistical descriptions of the data in log space. If you assume that 
log(N) is normally distributed about its average value, cycle life, N, will have a lognormal 
distribution. 


For our example, we will arbitrarily select 145 ksi for a stress level and calculate the average 
cycle life and standard deviation for that stress. 

H y = a+b log( 145,000) = 7.5 193 . 

The standard deviation is considered to be constant regardless of stress for this example, 



(6) The last step is to fit a probability distribution to the data, such as a normal or 
Weibull distribution. For this example, assume a normal distribution using the 
mean and variance calculated above. This will result in a lognormal distribution of 
cycle life. 
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The procedure outlined contains the basic steps used by the PFA fatigue life software. The 
PFA program, MATGRM, performs a least-squares curve fit and then fits a Weibull distribution to 
cycle life. The program then performs a Monte Carlo simulation at a specified stress level. 

To verify the PFA program, the six steps listed were followed. A lognormal distribution was 
assumed and 1,000 Monte Carlo simulations performed. The distribution of cycle lives was then 
plotted as a histogram. The same Waspaloy data were input into MATGRM with 1,000 simulations 
and plotted. Results are compared in figure 8. 
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Figure 8. Probability density function of Waspaloy data using lognormal distribution 

and the JPL PFA Weibull distribution. 


The histogram comparison verifies that a least-squares curve fit with a lognormal distribution 
is similar to JPL’s program with a Weibull distribution. Differences between lognormal and Weibull 
distributions can be debated, but they both accomplish the intended job; to put a distribution on life 
with only positive values allowed. 


B. Strain/Life Statistical Analysis 

Plastic strain fatigue data can be treated the same as stress fatigue data, in some cases, by 
linearizing in log space with the power law relation. 



where. 


A = plastic strain life coefficient 
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b = plastic strain life exponent 
N = number of fatigue cycles 
Ae p = plastic strain range . 

Hence, strain fatigue data can be curve fitted with the least-squares method as previously 
described. A linear curve fit in log coordinates may be sufficient for the data. In fact, the Waspaloy 
data in table 2 and the matching curve fit in figure 9 show this to be the case. 

Table 2. Waspaloy strain life fatigue data. 


Cycles, N 

Percent Strain, Ae 

log (N) 

log(Ae) 

74 

3 

4.30 

-4.20 

186 

2.5 

5.23 

-4.38 

867 

2 

6.77 

-A. 61 

12 

4 

2.48 

-3.91 

6 

6 

1.79 

-3.51 

48 

3 

3.87 

-4.20 

44 

3 

3.78 

-4.20 

61 

3 

4.11 

^1.20 

40 

3 

3.69 

-4.20 

586 

2 

6.37 

—4.61 

348 

2 

5.85 

-4.61 

363 

2 

5.89 

—4.61 

390 

2 

5.97 

-4.61 

119 

3 

4.78 

-4.20 

145 

2.5 

4.98 

^1.38 

592 

2 

6.38 

-4.61 

1,396 

1.5 

7.24 

-4.89 

5,086 

1 

8.53 

-5.30 

20,000 

0.75 

9.90 

-5.59 


In some cases, however, the strain data are nonlinear when plotted in log coordinates. A 
nonlinear curve fit in log space is accomplished by separating the strain into elastic and plastic 
components and curve fitting each separately. The total strain relation then becomes, 

^ = ^ + ^ = ANb +CN d , 

where the curve fit constants A, b, C, and d are determined by performing least-squares curve fits to 
the elastic and plastic portions of strain. 
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Waspaloy Strain Life Data 



log(N) 

Figure 9. Curve fit of Waspaloy strain life data. 


C. Conclusion/Comments 

The PFA fatigue life program was evaluated on the basis of its understandability, ease of 
use, and application to a typical fatigue problem. 

Theory for the materials characterization program is found in section 2.1.2 of the PFA core 
document. 4 The theory section has several major omissions: (1) no discussion of the basic 
fundamentals of fatigue or its statistical analysis, and (2) no discussion of the lognormal distribution 
for cycle life. The lognormal distribution should be the default distribution, yet the PFA program does 
not even have it as an option nor any discussion as to why the lognormal was omitted. The advanced 
level of mathematical and statistical notation makes it hard to get to the core of what the program 
really does. 

The materials characterization program structure is described in two separate sections, one 
titled “Materials Characterization Software” and another titled “Materials Characterization Pro- 
gram.” There is some redundancy here because both are describing the same program. The first 
section, 4.1, consists of a detailed, 60 pages, flowchart, with a brief written description of the 
program and major subroutines. The other section, 7.3, gives another flowchart and description of 
variables. Documentation is acceptable, but hard to follow and understand. 

A user's guide is provided in section 6.3, “Materials Characterization User’s Guide.” This 
section is useful for actually running the program. It does a good job of describing the input and 
output with examples provided. A program listing is provided in section 7. 3. 1.4. The PFA materials 
characterization program consists of 6,707 lines of FORTRAN code and 41 subroutines. 
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V. RESPONSE SURFACE METHODS— ANALYSIS OF A CANTILEVER BEAM 


The increased use of probabilistic structural analysis methods has initiated a need for a cost- 
effective and efficient way of interfacing with the finite element method. The main reason being that it 
is not pragmatic to conduct thousands of finite element runs. Therefore, response surface methods 
have been proposed as the link between probabilistic methods and the finite element method. In 
addition, response surface methods have been available for many years and are commonly applied in 
the pharmaceutical, chemical, and testing industries to characterize their respective outputs. 
References 7 through 1 1 provide a good background on the application and theory of such methods. 


A. General Background 

The basic concept of a response surface method as applied to finite elements consists mainly 
of fitting a chosen equation, usually a linear or quadratic polynomial, to a relatively small number of 
planned finite element runs. The number of runs needed depends on the number of parameters 
assumed to affect the response and the number of levels through which each parameter varies. For 
finite element applications, all of the parameters considered are required to be quantitative and 
continuous within the region of interest. 

The purpose of the research being presented in this section was to investigate the utility and 
accuracy of applying response surface methods to finite element results. For this purpose, response 
surface analyses, using the finite element method, were performed on a cantilever beam model under 
elastic and elastic-plastic conditions. In addition, two response surface designs were chosen for 
these analyses, and the results were compared. A cantilever beam was selected because of its 
simplicity and availability of a closed-form solution in the elastic case, and readily derivable iterative 
solution in the plastic region. 

The first step in conducting a response surface analysis, once the problem has been defined, 
is to determine the number of parameters to be considered. For the case of the cantilever beam, 
three parameters were assumed, namely. Young’s modulus, beam tip load, and beam length. Table 3 
shows the minimum, nominal, and maximum values chosen for each one of the parameters. 

Table 3. Parameter values for cantilever beam. 



Beam Length ( L ) 
i (in) 

Tip Load ( P ) 
(lb) 

Young’s Modulus 
(E) (lb/in 2 ) 

Minimum 

9.9 

855 

2.70E+07 

Nominal 

10 

900 

3.00E+07 

Maximum 

10.1 

945 

3.30E+07 

Percent off nominal +/- 

10 

5 

10 


The second step in the process is to select the equation to be fitted and to design an experi- 
ment matrix appropriate for the equation chosen. There are several response surface design matrices 
available in the literature to fit first- and second-order equations. The two particular designs 
selected for the cantilever beam example were a 2 3 factorial design, 7 and a Box-Behnken design. 8 11 
The first design consists of three parameters, each varied at two levels, maximum and minimum, and 
used to fit the following equation: 
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F = k 0 +k l L+k 2 E+k 3 P , 


( 1 ) 


where F is the value of response, L is the length of the beam, E is Young’s modulus, P is the tip 
load, and k 0 , k x , k 2 , and *3 are the constant coefficients that result from the response surface 
analysis. There are a total of eight finite element runs needed to fit equation (1). Table 4 shows the 
combinations of the parameter values needed for each run. The second type of design chosen for the 
cantilever beam example, the Box-Behnken design, consists of varying parameters to three levels; 
minimum, nominal, and maximum. In addition, this kind of design is used to fit a quadratic polynomial 
equation, such as: 


F = to+t]L+t2E+tjP+t4EL+t$PLi+tf l EP+tjL %E +t^P , ( 2 ) 

here to to tg are the constant coefficients to be determined. The Box-Behnken design requires a total 
of 13 finite element runs. Table 5 shows the value combinations for the cantilever beam using this 
kind of design. 


Table 4. Cantilever beam parameter values for 2 3 design. 


Runs 

Beam Length (L) 
(in) 

Young’s Modulus 
(E) (lb/in 2 ) 

Tip Load ( P ) 
(lb) 

1 

9.9 

2.70E+07 

855 

2 

10.1 

2.70E+07 

855 

3 

9.9 

3.30E+07 

855 j 

4 

9.9 ~i 

2.70E+07 

945 

5 

10.1 

3.30E+07 

855 

6 

10.1 

2.70E+07 

945 

7 

9.9 

3.30E+07 

945 

8 

10.1 

3.30E+07 

945 i 


Table 5. Cantilever beam parameter values for Box-Behnken design. 


Runs 

Beam Length (L) 
(in) 

Young’s Modulus 
(£) (lb/in 2 ) 

Tip Load (P) 
(lb) 

1 : 

9.9 j 

2.70E+07 

900 

2 

9.9 

3.30E+07 

900 

3 

10.1 

2.70E+07 

900 

4 

10.1 

3.30E+07 

900 

5 

9.9 

3.00E+07 

855 

6 

9.9 

3.00E+07 

945 

7 

10.1 

3.00E+07 

855 

8 

10.1 1 

3.00E+07 

945 

9 

10 

2.70E+07 

855 

10 

10 

2.70E+07 

945 

11 

10 

3.30E+07 

855 

12 

10 

3.30E+07 

945 

13 

10 

3.00E+07 

900 
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For either of the designs discussed above, the output to be obtained from the finite element 
runs is F ( F could be displacements, stresses, strains, etc). Also, note that the number of finite 
element runs required in either design is greater than the number of unknowns (K’s) to be 
determined. Therefore, a method such as least squares is used in order to solve for the constant 
coefficients. Once the least-squares approximations are complete, the response surface equation is 
fitted and one can proceed to verify its accuracy. In the case of the cantilever beam example, the 
fitted surface results of both designs were compared to theoretical results. 


B. Elastic Case Study 

The ANSYS finite element code Rev. 4.4 A (ref. 12) was used to perform the runs. A model 
using the three-dimensional (3-D) 20-node isoparametric solid element type consisting of 80 
elements and 621 nodes was built. This model is shown in figure 10 and the critical locations are 
identified for stress and displacement. The relatively fine mesh for this problem was selected to 
prevent possible model inaccuracies due to mesh density, thus eliminating mesh density as a 
parameter. 



Figure 10. Cantilever beam finite element model for elastic solution. 

The design matrices, tables 4 and 5, for each design were used to perform the finite element 
runs. The least-squares approximations of the constant coefficients for the designs were also 
obtained. The values of these coefficients are shown in table 6 for stresses ( s ), and displacements 
(y) at the critical locations. 

In order to check for the accuracy of the fitted polynomials with the coefficients of table 6, 
closed form solutions were obtained from reference 13 for the displacements and stresses at the 
critical locations. The beam tip displacement theoretical solution is 

y = PL 3 /OEI) , (3) 

where, 

I=( l /l2)bh 3 . (3a) 

I is the moment of inertia, b the base of the beam, and h is the height. From figure 10, b = h = 1, 
therefore, /= i /n. The maximum stress from the theoretical solution is 

5 = Mc/1 , (4) 
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Table 6. Constant coefficients for elastic analysis. 


Coefficient 

Linear Fit (Equation (1)) 

Quadratic Fit (Equation (2)) 

y-Disp (Node 186) 
(in) 

5-Stress (Node 2) 
(lb/in 2 ) 

y-Disp (Node 186) 
(in) 

s-Stress (Node 2) 
(lb/in 2 ) 

0 

2.423E-01 

-5.306E+04 

2.676E-01 

-8.953E+01 

1 

-3.644E-02 

5.306E+03 

-4.284E-02 

2.689E+01 

2 

4.072E-09 

2.133E-18 

-4.057E-09 

1.092E-07 

3 

-1.357E-04 

6.050E+01 

1.507E-04 

1.437E+00 

4 



1.215E-09 

-1.325E-13 

5 



-4.208E-05 

5.905E+00 

6 



4.523E-12 

-6.487E-17 

7 



3.978E-04 

-1.722E+00 

8 



-1.348E-16 

-1.820E-15 

9 



-4.23 IE- 10 

1.129E-05 


where. 


M = PL . 


(4a) 


M is the applied moment and c is the distance from the neutral axis to the extreme fiber. From figure 
10, c = 1 / 2 . By making all the substitutions into equations (3) and (4) then. 


and 


y = APL? IE , 

(5a) 

s = 6PL . 

(5b) 


Equations (5a) and (5b) are written in terms of the three parameters chosen for the response 
surface analysis. In addition, the percent error between the theoretical and the fitted surfaces was 
computed as follows: 


and 


y percent error = 100*(y tab , e 6 ->'eq( 5 o)) / yeq( 5 a) > 
5 percent error = 100%v table6 -5 cq(5b) )/5 eq(5fc) • 


(6a) 

(6b) 


Figure 1 la displays the y percent error surface (equation (6a)) using the linear fit coefficients 
(equation (1)) from table 6 for a constant Young’s modulus of 33e6 lb/in 2 . This surface provides the 
greatest error obtained for all moduli (£) between the range of 27e6 and 33e6 lb/in 2 . Therefore, only 
the modulus for 33e6 lb/in 2 is shown here. In addition, this surface displays the y percent error 
dependent upon the variation of beam length L and beam tip load P within the ranges chosen in figure 
10. The maximum error obtained using this linear fit for the tip displacements was about 1.5 percent 
and occurs on a beam 10. 1-in long, carrying a load of 945 lb. Figure lib shows the same information 
as in figure 1 la with the difference that in figure 1 lb the coefficients from table 6 for the quadratic fit 
of equation (2) were used instead of the coefficients for the linear fit. Figure lib indicates a 
maximum error of about 0.82 percent occurring on a beam 9.9-in long under a load of 855 lb. 
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For the case of the s percent error, figures 12a and 12 b display the errors associated with the 
constant coefficients for 5 from table 10 for the linear and quadratic fits, respectively. It is important 
to note that the dependence on Young’s modulus is zero. This is obvious in the stress equation {5b) 
and was also confirmed in the results of the response surfaces. Therefore, figures 12a and 12 b are 
valid for all values of Young’s modulus. For the case of the linear fit, shown in figure 12a, the 
maximum error in s is about 0.90 percent and occurs on a beam 9.0-in long with a load of 945 lb. The 
surface for the quadratic fit (fig. 12 b) indicates a maximum error of about 0.86 percent on a beam 




9.95 


Figure 12a. S percent error from linear fit response versus theoretical results. 
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It can be concluded that in this elastic analysis both response surfaces adequately fit the 
cantilever beam problem having length, Young's modulus, and tip load as parameters with the 
variations indicated in table 3. In addition, it is further concluded that the difference between the 
linear and quadratic fits for this particular example were negligible. This would make the linear fit 
more desirable since it requires only eight runs. Also, the coefficients for the linear fit (equation (1)) 
provide the sensitivity of the response due to each parameter because they are the partial 
derivatives dF/dL, dF/dE, and dF/dP. However, before a general recommendation can be made on 
specific response surface designs for elastic analyses, more work needs to be done. For example, 
research is necessary on the maximum range of the parameters for which the response surfaces 
would continue to be satisfactory for different kinds of problems. Also, response surface analyses 
need to be performed on problems with high stress concentration that have available closed-form 
solutions. These two areas, among others, would help in providing a better indication of the 
adequacy of specific response surface designs for general finite element modeling use. Intuitively, it 
seems that no one design would serve for all problems unless the explicit form of the relationship 
(i.e., inverse, exponential, etc) between the parameters is known a priori. 


C. Elastic-Plastic Case 

The ANSYS finite element code Rev. 4.4a was also used to perform the elastic-plastic runs. 
For the model shown in figure 10, tensile yielding begins at node 2 or any node at the same end and 
at the outermost fibers. However, ANSYS only recognizes yielding when the equivalent stresses at 
the integration points of the elements are greater than the yield stress. Since, the integration points 
are located a distance away from the corner nodes of the elements which contain node 2, the 
outermost fiber yields before ANSYS establishes yielding at the integration points. Therefore, the 
mesh of the model shown in figure 10 had to be refined with an increased density of elements at the 
extreme fibers. The refined model is shown in figure 13. In addition, the refined model is composed of 
552 eight-noded 3-D isoparametric solid elements and 936 nodes. The response designs used for 
the elastic analysis were retained for the elastic-plastic analysis and the values shown in tables 3 
through 5 were applied to the refined model. Plasticity was introduced by setting a yield stress, s y> of 
50,700 lb/in 2 and by choosing the bilinear kinematic hardening option in ANSYS with the stress- 
strain curve shown in figure 14. Finally, the least-squares approximation of the constant coefficients 
for the displacements, y, stresses, s , and, total strains, et, at the beam’s critical locations and for the 
first- and second-order designs are shown in table 7. 



Figure 13. Cantilever beam finite element model for plastic analysis. 
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stress 



Figure 14. Stress/strain curve for plastic analysis. 

Table 7. Constant coefficients for finite elements of elastic-plastic case. 


Coeffi- 

cient 

Linear Fit (Equation (1)) 

| Quadratic Fit (Equation (2)) | 

y-Disp 
(Node 186) 
(in) 

s-Stress 
(Node 2) 
(lb/in 2 ) 

ct-Strain 
(Node 2) 
(in/in) 

y-Disp 
(Node 186) 
(in) 

A-Stress 
(Node 2) 
(lb/in 2 ) 

cf-Strain 
(Node 2) 
(in/in) 




50700/E+ 



50700/E+ 

0 

2.451E-01 

2.014E+04 

-1.027E-03 

2.138E-01 

4.292E+05 

7.657E-03 

1 

-3.660E-02 

1.549E+03 

5.066E-05 

-3.741E-02 


2.528E-04 

2 

4.071E-09 

-2.292E-05 

1.744E-13 

-4.040E-09 

M W53EM 

1.127E-10 

3 

-1.372E-04 

1.843E+01 

6.029E-07 

2.1 19E-04 

-5.843E+02 

-2.497E-05 

4 




1.214E-09 

-4.387E-05 

-9.972E-12 

5 




-4.517E-05 

3.397E+01 

1.432E-06 

6 




4.516E-12 

-5.083E-07 

3.947E-15 

7 




2.612E-04 

3.960E+01 

-5.859E-05 

8 




-1.349E-16 

7.068E-13 

-2.653E-19 

9 




-1.783E-08 

1.546E-01 

6. 197E-09 


To verity the accuracy of the lilted response of the elastic-plastic beam, a theoretical solution 
was derived by the using the concepts described by Phillips 14 and Johnson and Mellor. 15 Figure 15 
shows a rectangular bar under symmetrical loads P . This figure can be idealized as two cantilever 
beams where one end of each beam represents the rigid connection of the other. The solution derived 
tor this case using a linear hardening stress/strain curve such as figure 14 is. 


and 

for small y 


where 


2E t E *bh 3 k 3 +( 3£ 2 (E-E,)s^)h 2 -3£ 3 P(L-x) )k 2 +( (E-E,)s 3 b) = 0 , 

k = \/R = (d 2 yUlx 2 )/( 1 +(dy/dx) 2 ) m = d 2 y/dx 2 , 
e t = yk , 

a- = Ee e +E t {e ,-e e ) , 
e e -s y !E . 


(la) 

(lb) 

(7c) 

(Id) 
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Figure 15. Rectangular bar under load P (symmetrical about o ). 


This solution was iteratively solved for each of the input combinations shown in tables 4 and 
5, and linear and quadratic response surfaces were fit for displacements, stresses, and total strains. 
However, by using this procedure, the verification in the plastic case became the comparison of two 
response surfaces; one for the finite element results and the other for the theoretically derived 
iterative results. This approach undermined the evaluation of the response surface method for the 
elastic-plastic case. Nevertheless, the constant coefficients for displacements, stresses, and total 
strains for the theoretically derived solution are shown in table 8. 

Table 8. Constant coefficients — theoretical solution of elastic-plastic case. 



Linear Fit (Equation (1)) 

Quadratic Fit (Equation (2)) 

Coefficient 

y-Disp 
(Node 186) 
(in) ' 

^-Stress 
(Node 2) 
(lb/in 2 ) 

ef-Strain 
(Node 2) 
(in/in) 

y-Disp 
(Node 186) 
(in) 

5-Stress 
(Node 2) 
(lb/in 2 ) 

^r-Strain 
(Node 2) 
(in/in) 

0 

2.454E-01 

-4.829E+04 

-2.804E-04 

2.083E-01 

-5.053E+05 

3.539E-03 

1 

-3.655E-02 

6.417E+03 

1.857E-04 

-3.654E-02 

-3.444E+03 

7.859E-06 

2 

4.041 E-09 

— 1.31 6E-05 

-5.696E-1 1 

4.072E-09 

-5.238E-05 

-5.744E-1 1 

3 

-1.361E-01 

3.889E+01 

2.054E-06 

2.167E-04 

1.197E+03 

-4.474E-06 

4 




1.212E-09 

-8.448E-05 

-5.639E-12 

5 




-4.512E-05 

-9.873E+01 

5.351E-07 

6 




4.483E-12 

-9.358E-07 

-6.239E-14 

7 




2.221E-04 

■ftMisniMi 


8 




-1.340E-16 

2.897E-1 1 

1.884E-18 

9 




-1.966E-08 

-8.222E-02 

1.690E-09 


The comparison of the results of the two response surfaces are presented here only for 
completeness. The percent error between the two surfaces (finite elements versus theoretical) were 
computed in a similar fashion as equation (6). 
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y percent error = 100*(y tab , e 7 -y table stable 8 , 
*t percent error = 100*(c„ able 7 -* ftable gV*i tables , 
s percent error = 100*(s table 7 -s table 8 )/s table 8 . 


(8a) 


m 

(8c) 

Figures 16a and 16b show the y percent error of equation (8a) for the linear and quadratic fits, 
respectively. The y percent error for either curve is relatively low and the difference between the fits, 
as in the elastic case, is negligible. In addition, figures 17a and 17b show the e, percent error of 
equation (8b) for both fits. Again, the difference between fits is negligible. However, the percent 
error between the finite elements and the theoretical solution indicates a maximum of about 5- 
percent error for figure 17a and 6 percent for figure 17b. This error could be a stack up of errors 
caused by the creation of both response surfaces plus the inherent finite-element modeling error. 
Finally, the s percent error of equation (8c) is shown in Figures 18a and 18b, where similar results, 
as in the total strains, are indicated. The only conclusion that can be drawn from the elastic-plastic 
analysis regarding the response surface method is that the differences between the linear and 
quadratic fits were negligible. 



Figure 16a. y percent error from linear fit response versus theoretical plastic results. 
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Figure 16 b. y percent error from quadratic fit response versus theoretical plastic results. 



Figure 17a. e, percent error from linear fit response versus theoretical plastic results. 




Figure Mb. e, percent error from quadratic fit response versus theoretical plastic results. 
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Figure 18 b. S percent error from quadratic fit response versus theoretical plastic results. 


VI. HPOTP DEMONSTRATION EXAMPLE 
A. Introduction 

This effort was performed as an example of incorporating a large finite element model into a 
probabilistic analysis. An ANSYS structural model is used to predict material strain at a stress 
concentration on a Waspaloy turbine disk. Model strain is compared with strain life data to estimate 
the number of mission duty cycles to failure. Methods generally described as response surface 
analysis (multiple regression) are used to formulate a closed-form equation for material strain as a 
function of temperature and yield strength. Least-squares methods are employed to curve fit cyclic 
strain-life data. A Monte Carlo simulation is performed to predict the distribution of mission duty 
cycles to failure. After a detailed description of the individual steps, a concluding summary is given. 

1. Problem . Six HPOTP first-stage turbine disks were found with cracks at the interstage 
seal pilot rib. 6 A sketch of the disk and location of cracking is shown on figure 19. The HPOTP disk 
undergoes a severe thermal shock during a normal mission duty cycle. When the pump shuts down, 
hot gas (456 to 940 °F) is replaced with hydrogen gas; within a few seconds, the gas temperature 
drops to about 0 °F. A typical thermal profile is shown in figure 20. The large drop in gas temperature 
causes a significant thermal gradient in the Waspaloy disk, resulting in a high strain gradient due to 
the coefficient of thermal expansion. As shown in figure 19, there are small (0.030-in) radii on the 
seal pilot rib. Cracks naturally originate in the radii due to a strain concentration effect and low cycle 
fatigue (LCF). For this problem, failure is defined as the initiation of a crack. 

2. Description of Analysis . The purpose of this analysis, whether deterministic or prob- 
abilistic, is to predict the number of mission duty cycles to failure. This analysis involves several 
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Note: Pilot rib undercuts required for access to grind curvic teeth 
Figure 19. HPOTP first-stage turbine disk. 



Figure 20. Hot gas temperature-oxidizer preburner. 

steps: (1) A thermal model is run to generate time-dependent temperature profiles for the disk. These 
temperature profiles are different, depending on the assumptions made for engine speed and gas flow. 
(2) Temperature profiles generated from the thermal model are then used as input to a structural model. 
The structural model also uses material properties and includes the effects of geometry. Material strains 
due to the thermal gradient are output. The structural model is a large, solid element model with 
nonlinear, plastic strain capability. A wire frame plot is shown in figure 21. (3) Strain data from the 
structural model are then used with strain-life cyclic test data to estimate the number of cycles to failure. 
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Figure 21. ANSYS 15° structural model. 


B. Probabilistic Approach 

As just described, this problem could become very complex if uncertainties in temperature, 
strain, material properties, and modeling error are fully considered. In order to work the problem in a 
reasonable amount of time, the number of input variables which were allowed to vary as inputs to the 
finite element model were reduced to temperature and yield strength. These variables were perturbed in 
a somewhat arbitrary manner, the purpose being to demonstrate how a response surface might be 
generated. Due to the model size, only four perturbation runs were made. 

Another simplification is that of distribution assumptions for temperature and yield strength. 
Both random variables are assumed to have the lognormal distribution due to the multiplicative form of 
the response surface and the need to restrict strain to values greater than zero. 

Uncertainties in cyclic life are accounted for with a least-squares fit to strain-life data which is 
linear in log space. A nonlinear least-squares fit is possible. More information on this subject can be 
found in references 16, 17, and 18. 

1. Response Surface A response surface was fit to strain data from four model runs. Strain data 
with corresponding temperature and material strength are given in table 9. The information was 
generated by Swanson and Goode. 19 

Table 9. Disk strain versus temperature and yield strength. 



Output Strain 

Input Temperature 

Input Yield 

Run No. 

£r(%) 

log £t (%) 

Temperature 

(R) 

log 

Temperature 

Sy 

(ksi) 

log Sy 

1 

1.605 

0.4731 

916 

6.82 

140 

4.942 

2 

3.169 

1.1534 

1,400 

7.244 

140 

4.942 

3 

6.892 

1.9304 

1,400 

7.244 

125 

4.828 

4 

2.581 

0.9482 

1,400 

7.244 

150 

5.011 
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Two forms of response surfaces were fit to the data, one linear and one nonlinear. The linear fit 


is. 


e = Po+P \(T-T q )+P 2 (W , (9) 

where /Jo is mean value strain and pi and Pi are constants. Using least-squares multiple regression, we 
solve for the constants. Note there are three unknowns and four sets of data. Least squares allows use of 
all the available data. The calculated values of constants are, 

A) = 3.5618 , pi =0.004776 , 

The nonlinear response surface is, 

e=p 0 


Pi = -0.17842 


( 10 ) 


which can also be written as, 


log e = log p 0 +p , logT +P 2 logSy . (11) 

Again, using multiple regression, the constants are, 

po = 13.568, Pi = 2.031 , pz = -5A64 . 


Typically, the response surface utilized depends on (1) which equation best fits the actual model data, 
and also (2) the type of probability distributions selected for the random variables. Table 10 summarizes 
the two curve fits and shows that, overall, the nonlinear response surface is a better fit. 

Table 10. Response surface comparison. 


Temp. 

Mat Yield 

Linear Response Surface 

Model Percent 

Strain Strain Error 

Nonlinear Response Surface 

Strain Model Percent 

(Percent) Strain Error 

916 

140 

1.61 

1.605 

0.00 

1.605 

1.52 

0.06 

1,400 

140 

3.92 

3.169 

0.24 

3.169 

3.60 

0.12 

1,400 

125 

6.59 

6.892 

0.04 

6.892 

6.68 

0.03 

1,400 

Mean 

Value 

1,279 

145 

Mean 

Value 

138.75 

3.02 

2.581 

0.17 

Average 

error 

0.11 

2.581 

2.97 

0.13 

Average 

error 

0.08 


2. Statistical Properties of Random Variables . Statistical properties of the random variables 
temperature and material strength must be determined or estimated before a probabilistic analysis can 
be completed. The mean values used in this problem example are given in table 10. The numbers shown 
are the average values of four inputs. Although these values are not necessarily correct, they were 
utilized to work this problem. Normally, some data base or experience base is used to determine 
statistical properties. Similarly, the standard deviation values were assumed and are not based on any 
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real data base. The process of obtaining statistical data (mean and variance) is a fundamental concern for 
each random variable. Table 1 1 contains a summary of the statistical properties. Notice, COV, median 
value, and the other lognormal parameters are calculated from the mean and standard deviation. 


Table 11. Statistical parameters for random variables. 


Parameter 

Temperature (T) 

Yield Strength (Sy) 

Normal Parameters 



Mean, fi x 

1,279 

138.75 

Standard Deviation, a x 

200 

13.88 

COV, o x ln x 

0.1564 

0. 10 

Lognormal Parameters 



Median value, . = = 

Vl-tCOV, 2 

1,264 

138.06 

y = log[Med(*)] 

7.142 

4.928 

<Ty = logfl+COV/] 

0.02417 

0.00995 


Note: If y = log x is normal (// and a 2 ), then x is lognormal. 


3. Probability Distributions . The lognormal distribution was chosen for both independent 
random variables in this problem, since it models only positive values of the dependent variable, strain. 
This assumption is important, because the normal distribution model allows for negative strain values to 
occur during a Monte Carlo simulation. A distribution limited to positive values was necessary, since 
negative strains cannot be used with cyclic life data. To check the effect of choosing a different 
distribution on input random variables, normal and lognormal distributions were compared for 
temperature and yield strength (fig. 22). Notice, for random variables with positive values, the PDF’s are 
very similar. Model difference is more noticeable for the temperature variable, which has a greater 
coefficient of variation. 
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The problem of determining statistical parameters for strain is greatly simplified if all variables 
are assumed to have the same probability distribution. By assuming that all variables are lognormal, the 
mean and standard deviation for log(strain) can be determined in closed-form from statistical properties 
of the individual random variables. Rewrite equation (1 1) as, 

z = log £ = log/J 0 +/3 , logT +p 2 log5y. (12) 

Lognormal parameters for z are defined as, 

/i ; =io g (/j 0 f' , ’sy’’] , (i3) 

<7, = log[(l+COvy i (l4COV s *)'’ 5 ] , (14) 

where T and Sy are median values, and COV is the coefficient of variation. Equations (13) and (14) were 
verified by performing a Monte Carlo simulation of equation (12) with lognormal distributions on 
temperature and yield strength. Results are summarized in table 12 for 100 to 10,000 simulations. Note, 
the mean and variance stabilize and converge to the closed-form calculations. 

Table 12. Lognormal parameter comparison. 




Monte Carlo Simulations 

Parameters 

Equations (13) and (14) 

100 

1,000 

5,000 

10,000 

l*z 

1.1486 

1.1400 

1.1728 

1.1481 

1.1406 

° 2 z 

0.3968 

0.3957 

0.4103 

0.4093 

0.4034 


A lognormal statistical distribution for strain has now been established and verified. The 
lognormal PDF is, 


fxix) ~JTZ < 7 , x cxp 


(log x-^y 


2a\ 


(15) 


where x is the strain value. A plot of the lognormal strain distribution is shown in figure 23 and 
compared with a normal strain distribution for illustrative purposes. The normal statistical parameters for 
strain are calculated similarly to the preceding example. Note that with a normal distribution, the left- 
hand tail contains negative strain values and is truncated in figure 23. Behavior of a distribution in the 
tail regions is very important when trying to determine probability of failure. The lognormal distribution 
has a thicker right-hand tail and thus a greater number of high strain values. 


Note that the strain distribution given in figure 23 is the total strain predicted from a nonlinear 
finite element model; it does not necessarily cycle about a zero strain value. Strain values recorded from 
cyclic test data are generated from a rotating test specimen which do, in fact, cycle about zero strain, i.e., 
the ratio = -L The effect of strain ratio should be considered when comparing disk model 

strain with cyclic test strains. 
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Figure 23. Strain distribution lognormal versus normal. 


C. Strain-Life Statistical Analysis 

Strain-life analysis is similar to stress-life analysis. Section IV. A of this report describes 
statistical treatment of stress-life data. The cyclic strain-life data for this problem were treated with an 
equation of the form, 

NsA &f ' (16) 


where, 

A = plastic strain-life coefficient 
b = plastic strain-life exponent 
N = number of fatigue cycles 
Ae = total strain range. 

Strain-life data were fitted with the least-squares method. Table 13 lists the data as received from 
Rocketdyne. 

Having the strain-life data, the next step is to calculate least-squares parameters Sxx, Syy, and Sxy 
to obtain a straight line curve fit in log space. Equation (16) is written as, 

log N = log A+b l°g(^r) . 07) 

which has the form 

y = a+bx , (18) 
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Table 13. Waspaloy low-cycle fatigue test results. 


Specimen 

Number 

Number of 
Cycles (AO 

T-axis 

(ln(A0 

Strain Range, e * 
(Percent) 

X-axis 

(ln(e/2)) 

29.00 

74 

4.30 

3.00 

-4.20 

13.00 

186 

5.23 

2.50 

-4.38 

16.00 

867 

6.77 

2.00 

-4.61 

30.00 

12 

2.48 

4.00 

-3.91 

27.00 

6 

1.79 

6.00 

-3.51 

HC1 

48 

3.87 

3.00 

-4.20 

HC2 

44 

3.78 

3.00 

-4.20 

LC1 

61 

4.11 

3.00 

-4.20 

LC2 

40 

3.69 

3.00 

-4.20 

HC3 

586 

6.37 

2.00 

-4.61 

HC4 

348 

5.85 

2.00 

-4.61 

LC3 

363 

5.89 

2.00 

-4.61 

LC4 

390 

5.97 

2.00 

-4.61 

64.00 

119 

4.78 

3.00 

-4.20 

6.00 

145 

4.98 

2.50 

-4.38 

61.00 

592 

6.38 

2.00 

-4.61 

26.00 

1,396 

7.24 

1.50 

-4.89 

11.00 

5,086 

8.53 

1.00 

-5.30 

51.00 

20,000 

9.90 

0.75 

-5.59 


* Strain range is total strain range from compression to tension for strain ratio = -1. 


where. 


a = log N 

b = b 

^ io *(¥) • 

The following least- squares parameters are calculated from table 13 data. 

Sx x = 3.91 , 

Syy = 70.87 , 

Sxy = -16.33 , 


mean value slope, 
mean value intercept, 


b = Sxy/Sxx = -4. 17 , 
a = Tavg-fe*Xavg = -13.26 . 
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The next step is to determine the basic statistical parameters for cycle life, lo g[N). Here, \og[N] is 
assumed to be normally distributed. For a given strain level, e/2, the mean and variance of log [AT) are 
calculated: 


mean value cycle life, 


y = a+bx , 


variance. 


<7 2 = 


n- 2 


(19) 

( 20 ) 


For the strain life data in table 13, <7= 0.40. 

For this problem, it is assumed that the mean value slope, intercept, and variance are constant 
over the range of test data. We must be careful to remember the strain level used in equation (19) is e/2, 
not total strain. During a simulation, total strain from the response surface equation must be divided by 2 
before using it to calculate cycle life. 


D. Monte Carlo Simulation 

A Monte Carlo simulation is used to generate the distribution of cycle lives. Each simulation 
consists of the following steps, 

(1) Generate total strain value from lognormal distribution. 

(2) Divide total strain value by 2 and use equation (19) to calculate the mean value of log [N], 
y = IogC/V). 

(3) Use mean value log(/V) and standard deviation from equation (20) to generate a cycle life 
from lognormal distribution. 

(4) Store the result (cycle life) and repeat steps 1 to 3 until the maximum desired number of 
simulations is reached. 

Figure 24 shows the results of 5,000 simulations. The cumulative distribution is used to 
determine the probability of reaching a particular cycle life. A cumulative distribution is very easy to 
calculate from a Monte Carlo simulation. Simulation results, such as cycle life, are sorted in order of 
increasing cycle life and then saved for plotting. After sorting, each succeeding cycle life corresponds to 
a particular percentage of the total number of simulations. If there are 5,000 simulations, the lowest 
cycle life corresponds to 1/5,000 of the total population; the 100th lowest cycle life corresponds to 
l/50th of the population. Table 14 gives some specific percentages with corresponding cycle life. 
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Figure 24. Cycle life — cumulative distribution. 

Table 14. Number of cycles to failure versus percent of total disk simulations. 


Number of Cycles to 
Failure 

Percent of Total Disk 
Simulations 

Percent Probability of Reaching a 
Particular Cycle Life Before 
Failure 

= 1- Previous Column 

56 

50 

50 

88 

60 

40 

140 

70 

30 

245 

80 

20 

547 

90 

10 

993 

95 

5 


E. Summary/Comments 

The primary goal for this effort was to demonstrate how a finite element model could be utilized 
to develop a response surface which is then incorporated into a Monte Carlo simulation of material strain 
and cycle life. This work should not be considered in any way complete; however, it does serve as a 
starting point for understanding the difficulties encountered in a typical application. The following 
comments are made relative to this problem: 
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(1) Development of a response surface from a finite element model is not a trivial matter. This 
problem was greatly simplified in order to work through the steps in a reasonable amount of 
time. 

(2) The structural model was not verified by test data and therefore modeling error could not be 
accounted for. 

(3) The ANSYS model included material nonlinearity. This is probably not the most efficient 
way to generate data for a response surface in which many runs are needed. A better 
approach would be to assume linear behavior since there is only local yielding. Neuber's rule 
could then be used to estimate notch strain. 

4) The average temperature and strain values are probably not correct, but were used to work 
through the problem. More effort is needed to characterize the input variables and level of 
uncertainty. 

(5) There is the question of whether or not actual strains on the disk can be directly compared to 
test strain data. Disk strains are local strains in a .030-in radius and are thermally induced. 
Strain-life data are generated from a smooth, round fatigue specimen at ambient temperature 
and are produced mechanically. Another significant difference is the strain ratio. The strain 
ratio, £ m in/£max- for disk strains is not known, while the test strain ratio is -1. 


VII. ALTERNATE PROBABILISTIC APPROACHES 


A. Introduction 

The purpose of this section is to supply some of the basic theory on calculating failure 
probabilities. This is by no means an exhaustive treatment of the subject, rather it is intended to give the 
reader a general understanding of how probability values may be calculated and how these techniques 
may be applied to large structural models. 

Probabilistic structural analysis attempts to predict the reliability (probability of nonfailure) of a 
structure using known, or assumed, statistical distributions on structural parameters. In general, such an 
analysis requires defining a mode of failure (material yielding, excessive deflection, etc.) and 
constructing a model to predict the response by which failure is measured. For example, to find the 
probability of material yielding in a structure; define failure as occurring when applied stress exceeds the 
yield stress, then generate a structural model to calculate stress in the structure as a function of the 
distributed parameters (also called random variables). The model relates the probability distributions 
(i.e., probability density functions or cumulative density functions) of the random variables to the 
probability distribution of the model response. Knowing the probability distribution of the response 
variable, the probability of material yielding can be calculated. 

Although various methods for determining the probability distribution (and thereby the 
probability of failure, Pf) of the response exist, most of these methods become very complex as the 
number of random variables increases. Several of these methods will be illustrated herein by solving a 
simple example problem. The first method will be an “exact” approach involving the integration of the 
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probability density functions. Following this will be several approximate methods including simulation 
and some basic concepts used in the NESSUS code. 

The example problem to be used is that of a prismatic bar of square cross section in simple 
tension (fig. 25). 



Figure 25. Prismatic bar in simple tension. 

In this example, failure is defined as material yielding, i.e., when the stress due to p exceeds the 
yield stress of the material, or 


where 



s = yield strength 
p = applied force 

♦ 

b = cross section width and height. 

In order to keep the problem simple, p is chosen as a constant and b and s are random variables defined 
in table 15. 


Table 15. Distribution parameters for example problem. 


Parameter 

s 

b 

P 

type 

Weibull 

Normal 

constant 

P 

10 6 lb/in 2 

1.0 in 

800,000 lb 

a 

10 4 lb/in 2 

0.1 in 

0.0 


The left-hand side of the inequality is commonly called the “failure function” or the “limit state 
equation.” The failure function is generally designated as “g” and is written as 



where failure occurs when 


£<0 . 
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Consider for a moment the term as “applied stress/' 


then the failure function is 


g=s-Z a . 

Now assume that has a known PDF, since s and Z, a have the same units, both PDF’s may be plotted 
on the same axis (fig. 26). 



Figure 26. Probability density function of applied and resistive stress. 

The region where the two PDF’s overlap indicates the range of values for which strength is less 
than stress; also, the area of the overlapping region relates to the probability of failure. Since g is a 
function of random variables, g is also a random variable. If the PDF of g is known, then integrating the 
PDF of g for all values such that g < 0 yields P/ffig. 27). 



Figure 27. Probability density function of random variable g. 

The area of the shaded region is the probability that g < 0 (also denoted as P[g < 0] ), or that the 
material has yielded. Calculation of the shaded region in figure 27 is given by 
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0 


P[g< 0] = j f g (g)dg » 


where /^(g) is the PDF of g,and -«» to 0 is the failure region. Problems that have only one random 
variable with a known PDF, such as g, present few difficulties. 


B. Direct Integration 


Generally, the PDF of g (or of the single response variable) is not known. This is the case for the 

P_ 

b : 


example problem shown in figure 25, where the distribution of b, rather than is known. Also, the 


example problem contains two random variables, adding to its complexity (this is sometimes called a 
two-dimensional problem). One classical approach to this type of problem is the direct integration of the 
joint PDF over the failure region. A joint PDF is a multivariable function that gives a probability density 
value for specified values of the random variables. For n independent random variables, the joint PDF is 
given by 




X = • • • Jn) 


where ffai) is the PDF of each random variable X ) (known as the marginal distribution). Once the joint 
PDF is known, the probability of failure may be calculated from 



(Q is the failure region) 


h n 

i 


n fi(Xi)dx n . 

I = 1 


. dx 2 jdx l , 


where the limits of integration are defined by the failure function g. Notice that this integral is the multi- 
dimensional equivalent to the integral for P[^<0], 


The example problem has a two-dimensional joint PDF, which may be written as follows 


f%(x) — f(s,B)( s ’b) - f s (s) ■ f/j(b) , 


where the two marginal distributions are defined as 




s-y 


P~\ 


exp 




(/J = 2 , y= 980,869.4 , rj = 21,586.6 => a s = 10,000 


and 




a B </ln 


exp 


l-^ii 


A plot of the joint PDF for the example problem is shown in figure 28. 


(Weibull Distribution), 

Us = 1 , 000 , 000 ), 

(Normal Distribution) . 
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Figure 28. Joint probability density function for example problem. 

Figure 29 is a contour plot of figure 28, which shows how the g function divides the variable 
space into safe and failure regions (not to scale). Thus, calculating the probability of failure by direct 
integration for problems with two random variables requires finding the volume under the joint PDF for 
all values of the random variables where g < 0. (The procedure is the same for problems of more than 
two random variables, but it cannot be visualized in this way.) 

6 

1.04 10 

6 

1.03 10 

6 

1.02 10 

6 

1.01 10 

6 

1 . 10 
990000 
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Figure 29. Contour plot of joint PDF of example problem. 


Since s and b are independent, the integral over f2 may be written as 


b 2 $2 

p f =j J m-mdsdb 


J *2 

f B {b) j f s (s)dsdb 




The limits of integration are obtained from the limit state and from the basic restrictions on the 
individual (or marginal) PDF’s. Since /s(s) is a three-parameter Weibull distribution, which is defined as 
zero for all s <y, the lower limit is 


s — y , 


The maximum value which S may have while remaining within the failure region Q is obtained 
from the g function (used here as a limit function or “limit state”). As can be seen from the contour plot, 
as long as g < 0, both variables are in the failure region, therefore 



which leads to the upper limit on yield strength 


S 2 



Although/flf/?) is a normal PDF, where < b < the fact that s 2 must be >/places a restriction 
on the range of b. The reduced range of b is found by replacing s 2 with gamma and solving for b, which 
yields 



Including limits, the integral for Pf is 

s fWr p/b 2 

P f = j Mb) j fs(s)dsdb , 
-JpTy r 


and performing the integration over s yields 


JpFy 

P/= f Mb)[l-e~ z2 ]db , 
ifWr 
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where 



This last integral was solved numerically since no closed-form solution exists for the CDF of a normally 
distributed variable. The resulting probability of failure was calculated to be 


P f = 0.145861 . 

Although this is an “exact” technique, (the only error in the answer comes from the numerical 
integration) the effort required for problems with many variables is so great that other procedures are 
much more attractive. 


C. Fast Probability Integration (FPI) 

The techniques known as “fast probability integration” avoid the difficulties of direct integration 
by approximating a problem as one which has a simple solution. To illustrate this, first consider figure 
30. This figure shows a contour plot of a joint PDF. The marginal distributions for this joint PDF are 
known as standard normal distributions, which means they have mean values of zero and standard 
deviations of one. The coordinate system in which the joint PDF is plotted is generally referred to as 
standard normal space, or u space. This special PDF has a particular property which makes it useful. 
Specifically, integration over any region bounded by a linear function can be reduced to a single integral. 
This property holds true for any number of coordinates/dimensions in the particular u space under 
consideration. In other words, a problem with 20 random variables, a multinormal PDF, and a linear 
limit state can be solved by performing one integration. This is in contrast to evaluating 20 nested 
integrals using the direct integration technique. 



Figure 30. Standard normal joint probability density function. 
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To evaluate the single integral, -<*> is used as the lower limit. The upper limit is the minimum 
distance from the origin to the limit state. (The point on the limit state closest to the origin is generally 
called the design point, or most probable point and is written as a* in standard normal space.) In figure 
30, the limit state is the level curve g = 0, and the minimum distance from the origin is the value (5. In 
terms of the coordinates u\ and ui, the function g may be written as mui+b-u\ (where m is the slope and 
b is the y intercept). Because of the symmetry of the PDF, we can write, 

„ mu 2 +b ~P 

^(M,)0(M 2 >y« 1 JH 2 = J <tf(z)dz = <t>(-fl) , 

— ©O 

where mu 2 +b is the equation for g = 0 solved for u\, and 0 is the standard normal CDF. For this special 
case, problems with many random variables can be solved with very little difficulty. 

In the real world, few problems have multinormal joint PDF’s and linear limit states. However, it 
is possible to approximate a real problem using this format. One way to do this is summarized as 
follows. 

• Approximate each nonnormal random variable as an “equivalent normal.” This requires 
finding means and standard deviations such that 

and ■ 

• Use the mean, and standard deviations of the equivalent normals to transform the limit state 
function into standard normal space (i.e., form g(a) by substituting x t = rx J w J +/i,' into gCs)). 

• Approximate the transformed limit state as a first order Taylor series function. This requires 
calculating partial derivatives of g. The approximation has the form 

Note that in the first step £* refers to the design point in the original variable space, i.e.. 


Once the first-order Taylor series function in standard normal space is known, a new design point and p 
can be found, allowing the probability of failure to be calculated using the single integral mentioned 
above. 

Generally, the design point found on the linear approximation to the limit state will be close to, 
but not equal to, the actual design point on the true limit state g = 0. Because the true design point should 
yield the best estimate of the probability of failure, an iterative procedure is used to converge to u*. The 
following paragraph describes a simple iteration procedure for finding &*. 

For simplicity, the iteration procedure will be described in standard normal, or u space. Although 
any point in u space may be used as the starting point, the mean values are usually chosen unless a better 



52 



initial estimate for the design point is available. This description begins at an arbitrary step (step m) 
where the vector (or point) has been found. The next step is to calculate the value of the failure 
function at & (m) , i.e., g(n (m) )- This yields a constant value, go. which defines the level curve containing 
the point a (m) - Calculating the gradient at & tm) and dividing by its magnitude produces a unit vector (£f< m ) 
in figure 31) perpendicular to g - go at the point ip m \ The unit vector describes the direction of the 
vector ^ ,n+ ri, which will describe the next approximation point. The vector is calculated as the 
sum of two vectors, both having the direction described by (£ m \ The first is found by projecting 
onto This results in a vector describing the minimum distance from the origin to the linear 
approximation of the level curve g = go- The second vector is added to account for the fact that may 
not be on the true limit state (i.e., g(M (m) ) * 0). Adding the two terms identifies a new approximation 
point, and a new linear approximation can then be calculated. Thus, the calculation proceeds as follows; 

go = g(ii (m) ) (find So) 


„(m) = grad (g(a (m) )) 
llgrad(g(ii (m) ))ll 


(calculate ^ m) ) 


ii to+D = [a (m) + .g (m) (calculate next point) . 

Ilgrad(g(ii tm, ))ll 

Using n ( - m+l '> as the new approximation point, the iteration continues by forming a new linear 
approximation for the curve g = g(iL (m+l) ), calculating ^ m+1) , and The procedure is then repeated 

until the estimate of u* converges. After converging to &*, the safety index is calculated as 

p= 1 1 tf* 1 1 , 

and the probability of failure is 

P f =®(-P) . 
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Effectively, for a two-dimensional problem, this procedure calculates a line tangent to a level 
curve of the surface g at each iteration step (although only one is shown on figure 31). It then finds the 
vector of minimum length from the origin to the tangent line. Adding a correction factor to this vector 
brings the point of the vector closer to the limit state, g = 0. The vector plus the correction factor 
describes a new point in u space that is used to calculate a new line tangent to the g function. This 
process is repeated until it converges to the point on g = 0 closest to the origin (commonly called the 
most probable point or MPP). 

One of the more basic FPI techniques, the Rackwitz-Feisler (R-F) algorithm, is very similar to 
the procedure just described. A flow chart of the R-F algorythm is given in figure 32 (reproduced from 
Madsen, Krenk, and Lind, reference 20, page 98). The equations in this flow chart are written in terms of 
the original variables, zi , and do not explicitly involve standard normal space. 



Figure 32. Flow chart for Rackwitz-Feisler algorithm. 
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The example problem was solved with the R-F procedure using a FORTRAN program and also 
by a hand calculation. The results of these calculations are as follows: 

hand calculation (three iterations) Pj = 0.1471 
computer calculation Pf = 0.1441 . 


Although the R-F algorithm is fast, it is not well suited for problems where the limit state is 
poorly approximated by a tangent hyperplane (the term hyperplane is used when the number of random 
variables is arbitrary). Also, certain types of distributions are not represented well by a normal 
distribution. A more advanced algorithm, known as the Wu-Wirsching (W-W) algorithm, 21 employs 
techniques which minimize errors due to both of these sources. 


The W-W algorithm reduces the error resulting from the linear approximation of the limit state in 
the R-F method by approximating the exact limit state as a quadratic Taylor series function (mixed terms 
are ignored): 


*»-*“*>*£ E 




(x r x .*y 


In the above equation, is the R-F design point, which is on the limit state, therefore 


g(x.*) = o . 


After calculating the derivatives for the quadratic approximation, the quadratic limit state is transformed 
into a linear limit state using one of several strategically chosen transformations. This reduces the 
computational effort required to find the safety index (p) while retaining the information contained in 
the quadratic approximation. It may be of interest to note that since the limit state is initially 
approximated as a second-order polynomial, the W-W algorithm may be classified as a second-order 
reliability method (SORM). This is in contrast to the R-F algorithm that approximates the limit state as a 
first-order polynomial, and is referred to as a first-order reliability method (FORM). 

Errors caused by an equivalent normal which poorly approximates the original distribution are 
reduced by using a “three- parameter equivalent normal.” The three-parameter equivalent normal is 
simply a two-parameter normal distribution, multiplied by the scale factor A. 

three-parameter equivalent normal = A<f)(x) . 

The three parameters (A,/r,< r) are found using a least squares approach. Beginning with the R-F values of 
H and crof each variable, the linearized limit state is written in the form 

g(jc) = a 0 +aiX,+T = 0 , 

where Y is the sum of all terms other than X t . This reduces the problem to two random variables, where Y 
is known (from the initial R-F analysis) and we wish to determine A, cr,, /i,- for X v The unknown values 
are found by minimizing the square of the error between the approximate and actual distributions 
modified by an appropriate weighting function. This is formulated as minimizing E, where 

oo 

E = f [A<t>(z)<f>(v)-F x (x)<p(v)] 2 dx , 
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subject to the constraint 


ci o+ci jX j+Y — 0 , 


where 


x-Hi J x-ti y 
z = -c~ and v = -^— 

w I v V 


Although the three unknowns may be found by applying an optimization routine to minimize E, it is 
more efficient to require that the area under the two functions is equal over the failure region. 


A = 


f F x {x)<j) (v)dx 
Jn 

f <b{z)<t>{v)dx 
Jo 


Substituting this last equation into the equation for E reduces computational effort by reducing the 
number of variables to be found by the optimization routine. 

The general procedure for this algorithm is similar to the R-F algorithm. Each step begins with 
an R-F solution for ft and the corresponding design point and statistical moments. An improved solution 
is then obtained by using the quadratic limit state approximation and the three parameter equivalent 
normals. As with the R-F algorithm, the procedure is repeated until a convergence criteria is met. 

The example problem was solved using a FORTRAN program of the W-W algorithm. Although 
this algorithm has been shown to be more accurate than the R-F procedure, 21 no difference between the 
results of the two different solutions was observed. This is probably because the true response function 
of the chosen example is almost linear near the design point. 

Finally, the probability of failure was calculated using Monte Carlo techniques. Using a modified 
version of one of the JPL programs, the probability of failure was calculated to be 

P f = 0.14556 , 

for 1,000,000 simulations. 


Although many of the algorithms for probability estimation are more advanced than the ones 
presented here, two concepts are central to most. Briefly, these concepts are: 

(1) Using some type of transformation from “jc” space (original variables) to space 
(standardized variables) 

(2) Use of an approximate response function. 

Most probabilistic algorithms will use these concepts in some form. For more information on the basic 
techniques, reference 21, and a paper by Chen and Lind 22 is recomended. For more advanced 
techniques, the paper by Thacker and Wu 23 and its references should supply a good starting point. Also, 
a paper by Khalessi, Lin, and Trent 24 along with its references may be of interest. 
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D. NESSUS (Probabilistic Finite Element Software) 


NESSUS (numerical evaluation of stochastic structures under stress) is an automated probabilis- 
tic structural analysis program that is capable of response, reliability, and risk (cost of failure) assess- 
ment . 25 Reliability and risk analysis are also available for structural systems using a fault tree to relate 
the bottom events. One of the unique features of NESSUS is the automatic generation of "response sur- 
faces" from implicit (finite element) models. This ability allows the user to perform probabilistic analy- 
sis on very complex components and/or systems defined implicitly by either a finite element or 
boundary element model. Without this feature a probabilistic analysis of an implicit model would 
require the user to manually alter the model to generate the required explicit response function. This 
would be very time consuming. Also, the advanced mean value (AMV) method employs iterative 
procedures that will re-analyze the finite element model until reaching a user defined convergence 
tolerance. NESSUS allows finite element random variables defined as loads, material properties, 
element properties, geometry, or some combination of these. This gives the user great flexibility in 
describing the analysis problem. 

The NESSUS software package is arranged into several modules. While most modules can be 
run as a “stand-alone” analysis, they are designed to be stacked together to form an automated analysis 
procedure. Three of the modules may be considered the primary functions of NESSUS. These modules 
are PFEM, FEM/BEM, and FPI. 

The FPI module does the actual probabilistic calculations. A variety of probabilistic calculation 
methods are available, and include FORM, SORM, FPI, and several sampling techniques. During a 
typical analysis, FPI will receive random variable definitions from PFEM and response function 
information from FEM or BEM. FPI will then perform the selected analysis and (conditionally) return 
information to the controling module about how the analysis should continue. 

The PFEM module controls probabilistic analysis of a finite element or boundary element model. 
In this module, the user defines random variables for finite element quantities and perturbations for 
generating the response surface. PFEM also controls various options for the mean value (MV), AMV, 
and adaptive importance sampling (AIS) probabilistic analysis methods. In general, PFEM controls the 
interaction between the FEM/BEM and FPI modules and the analysis techniques used. 

FEM/BEM’s are the finite element and boundary element modules, respectively. These modules 
allow the user to define complex implicit models (used to generate approximate response functions) in 
either a finite or boundary element format. These modules add much versatility to NESSUS since the 
user would otherwise be restricted to closed-form response equations. 

Other modules in NESSUS include: 

SIMFEM, which allows simulation to be applied directly to a finite element model 

RISK is used to compute component risk in terms of cost of failure 

PRE is a preprocessor which “uncorrelates” field variables given means, std., and a correlation 

model 

SRA coordinates PFEM with fault tree analysis for system reliability and risk assessment. 
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Although FPI techniques are considered very efficient, an FPI analysis using an implicit response 
function (like a finite-element model) may still require much computer time. To develop the explicit 
response function used by FPI techniques, it is necessary to solve the finite-element model several times 
for different values of the random variables. Because solving large finite-element models can be costly, 
efficient probabilistic analysis requires an algorithm that minimizes the number of model solutions 
required. One such algorithm, the AMV method, is available in the NESSUS software package. The 
AMV method was developed specifically for analysis of implicit response functions and uses very few 
function evaluations. 26 


The AMV procedure may be used for various types of analysis. It may be used to calculate either 
response values for a given probability (called PLEVELS analysis) or probability for a specified 
response value (called ZLEVELS analysis). Also, the NESSUS user has the option of using either a 
linear or an incomplete quadratic (mixed terms are ignored) as an approximate response function. For 
simplicity, the AMV method is described for a PLEVELS analysis using a linear response function. 


A brief summary of the AMV procedure is described in the following steps: 

( 1 ) Construct a linear approximation to Z (the true response) at the mean values 


Z(X)=Z(u)+t H 


X * 


= z, 


(2) Using FPI techniques, find zo and the MPP such that 


P[Z X < z 0 ] = requested probabilities . 


This constitutes the MV solution. 

(3) Re-analyze Z at the MPP found from the MV solution and recompute zo, the new zo is the 
AMV solution 

(4) If desired, form a new Z\ about the AMV MPP and return to step 2 until zo converges. 

For a more detailed description of this procedure, consider a continuous nonlinear response 
function 

Z = ZQ 0 , 

where X is a vector of n random variables 

X=(X u X 2 ,...,X b ) . 

Using a Taylor series evaluated at the mean values, Z may be written as 
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- Z=Z(u)+i #1 [XrUt + HQ Q 

1 — 1 O/L | I j^ii 

= a 0 +ta i X i +HQ Q 
= Z ,+//(£). 

Where i± is a vector of mean values and HQ 0 represents the higher order terms. The linear approxi- 
mation of Z at the mean value is Zi- By ignoring the higher order terms, the coefficients for Zi may be 
found using perturbation techniques requiring only n + 1 function evaluations. Once the coefficients ao to 
a„ are known, an analysis is done based on Zi - This is known as the mean value or MV analysis. 


of Z\. 


For the PLEVELS procedure, the MV analysis begins by using the mean and standard deviation 

n 

M Z, = a 0 + a i it* i * 


cr| = i . 

- 1 i=i 

These moments are used to estimate a range of probabilities (pj through ptf) and the corresponding 
(z<J through z") that span the user requested probabilities. In other words, find a range of zo values such 
that 


P[Z, < z^,] < requested probabilities < P[Z, < zg] • 

This step generates approximate points on the CDF curve of Z\. FPI routines are then used to improve 
the po estimate for each zo value and interpolation is used to approximate a zo value for the user 
requested probabilities. With the interpolated zo values, FPI is used once more to calculate more precise 
values of po . Also, the calculation produces specific values of the random variables corresponding to 
each po. These values are known as the design points, or the MPP’s for each po. This completes the MV 
analysis. Note that although FPI was run several times, the analysis is fast since Z\ is explicitly defined. 
The vast majority of the computation effort lies in the generation of Z\ from the finite element model. 

The AMV procedure calculates improved response values by calculating values for the function 
H(X). This is done by re-analyzing the response function (finite element model) at the MPP’s found in 
the MV analysis. When the actual response function is nonlinear, H(X) will be nonzero, and is given by 

HQ*) = Z(£*)-ZiC£*) , 

where x* is the MPP. Adding the correction factor H(x*) gives an improved approximation to the true 
response function in the vicinity of &*. This in turn allows for improved estimates of po and the 
corresponding MPP. 
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An iteration procedure, known as AMV+, is also available in the NESSUS code. This procedure 
improves the estimates obtained from AMV by constructing a new Z L at the MPP found by AMV. The 
new Z\ is then used in a MV analysis and again in AMV, repeating this cycle until zq converges. 27 

Although the MSFC team has spent limited time reviewing NESSUS, several advantages and 
disadvantages of the NESSUS probabilistic finite element code have been noted and are listed below: 

Advantages 


• Automatic response surface generation and design point search 

• Random variables defined using mean and standard deviation for many types of distributions 

• Changing the type of probabilistic analysis is easy (FPI, Monte Carlo, etc.) 

• Format of input data is relatively simple 

• FPI algorithms produce sensitivity values as well as probability values. 

Disadvantages 

• The NESSUS finite-element package is limited in the types of elements it supports and in that 
only one element type may be used for each model 

• FPI algorithms may not be very accurate for problems with highly nonlinear response functions 
(results were compared to Monte Carlo results using explicit response equations) 

• The operator must guess as to what perturbation step size will give the best results 

• Defining geometry parameters as random variables can be quite time consuming when a large 
number of nodes are involved. 


VIII. PROBABILISTIC BASED DESIGN CODES 


Uncertainties in the definition of loads and environments, material properties, geometric 
variables, manufacturing processes, engineering models, analysis tools, and all types of testing including 
verification and certification lead to uncertainties in space vehicle and structure design, and ultimately 
safety. Clearly, quantifying and understanding “problem uncertainties” and their influence on design 
variables develops a better engineered, designed, and safer structure. Two formats are available for 
characterizing design uncertainties: (1) deterministic/safety factor and (2) probabilistic/ reliability. The 
classical deterministic analysis approach accounts for design uncertainties in "lump sum" fashion by 
multiplying the maximum expected applied stress (a single value) by a factor of safety. Often, design 
verification is achieved by applying a worst case loading (e.g., a 3-sigma load condition multiplied by 
the safety factor) to the structure and testing to failure. In contrast, the probabilistic format attempts to 
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map each individual parameter uncertainty into a probability density function. A test constructed data 
base gives the best characterization of random variables. If test data are not available, then the engineer 
must make assumptions concerning the parameter distributions. Once the distributions are defined, 
transformation equations are used to combine the density functions into a cumulative distribution 
function of the design variable; for example, applied stress. In this case, the design parameter has an 
uncertainty that is quantified in terms of risk. 

To develop a NASA structural design code based solely on the probabilistic format without 
compromising the structural safety of hardware design is not achievable for the near future. This 
statement is especially true with respect to the currently available probabilistic engineering analysis tools 
and test verification programs. In general, the analytical tools that have been developed are difficult to 
understand and implement into a design procedure; and more importantly, the methods have not been 
test verified or universally accepted by the engineering community. Before a probabilistic based design 
code or program can be successful, design engineers must develop an experience and education base in 
this field. Most engineering schools do not offer probabilistic based design courses as part of the 
curriculum. 

True reliability must be demonstrated not simply estimated from an engineering analysis. Until 
failure and failure rate data bases are available from experience, probabilistic methods can best be 
utilized as a design tool to help identify the sensitivities of problem parameters. Furthermore, 
"demonstrated structural reliability" is virtually an impossible task due to the expense and small number 
of structures that NASA builds. However, it may be feasible to develop a more consistent structural 
design code that uses the probabilistic format in combination with the accepted safety factor approach to 
design. The civil engineering profession has successfully used a combined format in the development of 
the load and resistance factor design (LRFD) code for steel structures. 28 Developing this concept for 
application within NASA offers a natural extension to the current MSFC probabilistic analysis team 
work tasks (long term objective), and provides a practical area for future research. 


IX. REPORT SUMMARY AND RECOMMENDATIONS 


In March of 1991, MSFC formed a special task group to review the probabilistic engineering 
analysis methodology PFA developed by JPL. 2 4 The JPL team was led by Dr. Nick Moore as part of the 
Engine Certification Project funded through codes Q and M (project start, 1985). A number of important 
findings have resulted from the Marshall team investigation and are documented within the body of this 
report. Comprehensive background information and technical details have been provided for the 
interested reader. The sections below provide: (1) a summary of the findings and (2) the 
recommendations of the task team review. Note, this report presents the work of the Marshall 
probabilistic analysis team during 1991-1992. The MSFC team members are; Dr. John Townsend, 
Mario Reinfurth, Rene Ortega, Charlie Meyers, Jeff Peck, and Bob Weinstock. 


A. Summary of Task Group Findings 

(1) The PFA methodology is not a generic “tool” that can be easily applied to most engineering 
problems. In almost every case, the software is problem specific and must be restructured for 
a given application. 
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(2) The JPL reports documenting the probabilistic approach and methods are extremely difficult 
to penetrate and understand. 

(3) The Monte Carlo front end feature of the PFA analysis program is the most valuable and 
useful part of the JPL work. However, the complexities such as Burr curve fitting and 
Bayesian updating, appear to be unnecessary. 

(4) The JPL team has not identified specific details of the engineering design process and how 
the PFA methods apply. In particular, the technique by which an engineer integrates 
complex models, such as finite elements, into a probabilistic framework is very unclear. 

Based on our findings, the JPL team has not developed an engineering tool that can be easily 
understood and applied. As part of this investigation, we have looked at a number of other methods to 
perform probabilistic engineering analyses. Some of these methods, namely NESSUS and FPI (Fast 
Probability Integration) offer a more practical probabilistically based tool for the engineer. 

It is important to note that the task group considers the introduction of probabilistic methods into 
design practices a valuable and necessary objective. Incorporation of distributions/uncertainties and 
sensitivity analyses into current design practices would aid in: (1) more carefully considering the basis 
and ranges for estimates of loads and material life, (2) avoiding overly conservative designs, and (3) 
characterizing the risk of failure. 

While the MSFC task group does not believe that the PFA methodology offers a practical, stand- 
alone, engineering tool for probabilistic analysis, we recognize the “ground breaking” efforts of the JPL 
team. In particular, the JPL effort has succeeded in introducing and advancing the probabilistic approach 
within NASA, with particular emphasis on coupling deterministic analysis with statistical distributions 
of significant parameters. 


B. Recommendations 

The MSFC team has been active in probabilistic design analysis for about 18 months. We have 
examined the PFA methods, other design tools such as NESSUS and FPI, conducted/participated in 
training sessions and seminars, and worked jointly with university researchers. Based on our 
investigations, the following recommendations are presented. 

1. PFA Methodology . The PFA methodology should not be adopted as a primary method for 
probabilistic structural analysis at NASA. With the exception of the Monte Carlo simulation front end 
feature, the analysis tools as currently constructed have limited application. The JPL front end 
simulation software is a useful FORTRAN product and should be made available to interested parties. 

2. NESSUS . The NESSUS probabilistic finite element software package is a good engineering 
analysis sensitivity tool for use in component structural design. If this package were more user friendly 
with broadened application features, NESSUS could become a standard design tool for probabilistic 
analysis in NASA applications. In particular, the current NESSUS code should be modified and 
expanded to allow for linking standard finite element codes (such as NASTRAN, ANSYS, etc.) and to 
allow for usage of multiple finite elements. 
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3. Investigate LRFD . Ten years ago, the civil engineering profession introduced probabilistic 
methods into structural design practices. The goal was to produce designs with a more consistent level of 
safety against structural failure. Their methodology is termed LRFD. Because LRFD moves beyond a 
deterministic factor of safety approach, it may offer beneficial insights and applicable methods. It is 
recommended that the MSFC task team investigate the applicability of LRFD. 

4. Near-Term Integration of Probabilistic Design Methods . For the near-term, the goal should 
not be to replace current engineering design practices with probabilistic methods. Rather, the goal should 
be to supplement current safety factor deterministic approaches with probabilistic methods. The 
emphasis should be on gradual introduction of probabilistic methods (1) to characterize and determine 
the effects of loads and material property uncertainties and (2) to estimate the risk of particular failure 
modes, such as crack propagation. 

In addition, the research efforts on probabilistic methods should not focus on the development of 
a “system reliability number.” While system risk quantification is desirable, more realistic and 
achievable progress can be made by focusing on sensitivity analysis of life drivers at the 
component/failure mode level. The latter emphasis will be of relevancy to the practicing design engineer 
and ultimately lead to a more robust design. 

5. Formulation of a Small Agencv-Wide Committee . To avoid duplicative effort across NASA 
centers and to achieve a consensus on direction and specific work efforts, it is recommended that a small 
working group be established. Initially, the committee would be chaired by NASA Headquarters, 
comprised only of NASA personnel, and be limited to a size of approximately 10. The committee would 
aid in charting probabilistic design project direction and funding. 

6. Survey Existing Probabilistic Design Efforts in NASA . It is apparent that there are a number 
of other engineering groups active in probabilistic design methodology at NASA, with little 
communication among groups. We recommend that a survey be conducted to determine : (1) what work 
is being done; (2) where and by whom; (3) the objectives, goals, and duration of the work; and (4) the 
cost. 


7. Training . The education and training game plan should be to develop expertise within a small 
group of NASA engineers with the purpose of these individuals developing and presenting probabilistic 
design courses. It is both ineffective and costly to continually contract out training for large numbers of 
NASA personnel. 
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